Rotating a Vector to Find a Star's Angle

Building the rete for a paper astrolabe yesterday, I kept coming back to how the thing actually works: you don’t calculate a star’s position so much as turn the disc until it agrees with the sky. Rotation is the whole operation. Line up Vega’s pointer with the current hour, and the tympan underneath reads back your time. The machine computes by moving.

Which is exactly what CORDIC does. If you’ve ever wondered how a 1980s pocket calculator found sin(40°) with a chip that couldn’t multiply — an HP-15C, an 8087 coprocessor doing its trig — this is the trick. You start with a vector pointing along the x-axis, then rotate it toward your target angle in a fixed set of ever-shrinking steps: about 45°, then 26.57°, then 14.04°, each one the arctangent of a power of two. Overshoot? The next step turns back. Every rotation is just a shift and an add, because multiplying by 2⁻ⁱ is free in binary. When you run out of the angle, the vector’s coordinates are your cosine and sine.

import Foundation
let K = 0.6072529350088813  // gain from the shrinking rotations
func cordic(_ theta: Double, _ n: Int = 16) -> (Double, Double) {
    var x = K, y = 0.0, z = theta
    for i in 0..<n {
        let d = z < 0 ? -1.0 : 1.0
        let f = pow(2.0, Double(-i))
        (x, y) = (x - d*y*f, y + d*x*f)  // rotate by ±atan(2^-i)
        z -= d * atan(f)
    }
    return (x, y)  // (cos, sin)
}
let (c, s) = cordic(40.0 * .pi / 180)
print(String(format: "cos=%.6f sin=%.6f", c, s))  // cos=0.766049 sin=0.642783

The same walk in Perl, which reached its own peak around the same time these chips did:

my $K = 0.6072529350088813;
sub cordic {
    my ($theta, $n) = @_; $n //= 16;
    my ($x, $y, $z) = ($K, 0, $theta);
    for my $i (0 .. $n-1) {
        my $d = $z < 0 ? -1 : 1;
        my $f = 2 ** -$i;
        ($x, $y) = ($x - $d*$y*$f, $y + $d*$x*$f);
        $z -= $d * atan2($f, 1);
    }
    return ($x, $y);
}
printf "cos=%.6f sin=%.6f\n", cordic(40 * 3.14159265 / 180);

That K out front bugged me until I worked out why it’s there: each rotation stretches the vector a little, and since the set of steps is fixed, the total stretch is a constant — 1.6468 — so you just pre-divide by it once. The engraver of an astrolabe does something similar. The stereographic projection is baked into the plate before a single star is ever sighted; at the eyepiece you only rotate. Both machines front-load the hard geometry so the nightly work is turning a dial. Sixteen steps got me to five decimals against a star sighting off by three minutes — the algorithm was the honest half.