Search code examples
javascriptcoordinate-transformationastronomy

Bug with astronomical coordinate conversion


As part of a bigger project, I'm developing a JavaScript library that will help me generate ephemerides for pretty much whatever you want. The current piece I'm working on is coordinate conversion. As a tester project, I'm generating a list of degrees ON the ecliptic (that is, longitude 34, latitude 0) and their corresponding horizontal azimuths for a given time and place. The problem I'm getting is that in the middle, between ecliptic degree 155 and 334, all the azimuth jump into the 20000 range. What's going wrong?

The HTML script:

function main() {
    var i, e, h, op;
    var olon = -121.46888;
    var olat = 38.55556;
    var jd = julian(1971,3,21,6,38,0);
    var ecl = getNuta(jd).ecl;
    op = "Equinoctal Az by Ecl Degree<br><br>"
    for (i = 0; i < 360; i++) {
        e = eclEqu(i, 0, ecl);
        h = equHor(e.ra, e.dec, olat, lst(jd, olon));
        op += "Ecl: " + i + " = Hor: " + h.az.toFixed(3) + "<br>"
    }

    document.write( op);
}

And the conversion functions:

function eclEqu(lambda, beta, ecl) {
    // Converts ecliptic to equatorial
    // -> 3x float in degrees
    // <- {ra:, dec:} 2x float in hours, degrees

    var l, b, e, a, d, up, dn;
    l = lambda * co.dtr;
    b = beta * co.dtr;
    e = ecl * co.dtr;
    up = Math.sin(l) * Math.cos(e) - Math.tan(b) * Math.sin(e);
    dn = Math.cos(l);
    a = Math.atan2(up,dn);
    d = Math.asin(Math.sin(b) * Math.cos(e) 
        + Math.cos(b) * Math.sin(e) * Math.sin(l));
    return {ra: a * co.rth, dec: d * co.rtd};
}

function equHor(alpha, delta, phi, lst) {
    // Converts equatorial to horizontal
    // -> 4x float, in hrs, deg, deg (obs. lat), and hrs
    // <- {az:, alt:} 2x float in deg

    var  d, p, ha, A, h, up, dn;
    d = delta * co.dtr;
    p = phi * co.dtr;
    ha = (lst - alpha) * co.htr;
    up = Math.sin(ha);
    dn = Math.cos(ha) * Math.sin(p) - Math.tan(d) * Math.cos(p);
    A = norm(Math.atan2(up, dn), 360);
    h = Math.asin(Math.sin(p) * Math.sin(d) 
        + Math.cos(p) * Math.cos(d) * Math.cos(ha));
    return {az: A * co.rtd, alt: h * co.rtd}; 
}

norm is a simple function that normalizes a number into the range of 0 to x (360 in most cases above) by adding/subtracting multiple of x until in the range:

norm(n, x) = n - ( x * int(n/x))

I applied norm to every output, which helped, but I still get these crazy jumps:

Ecl: 153 = Hor: 2.537
Ecl: 154 = Hor: 0.557
Ecl: 155 = Hor: 105.103
Ecl: 156 = Hor: 103.214

and

Ecl: 333 = Hor: 289.018
Ecl: 334 = Hor: 287.038
Ecl: 335 = Hor: 178.622
Ecl: 336 = Hor: 176.734

I'm also concerned - it seems there should still be an even, smooth distribution (even if it's not linear) around the horizon, but it doesn't seem to be the case.


Solution

  • I think the problem is the statement

    A = norm(Math.atan2(up, dn), 360);
    

    in the function equHor(). The result of atan2 is an angle in radians. You wouldn't probably normalize it with 360.