Skip to content
← All posts
javascripttypescriptnumericsperformancenpm

Converting NREL's Solar Position Algorithm to JavaScript: A Journey in Precision

4 min read

Porting a reference C algorithm to JS accurate to ±0.0003°. The catastrophic cancellation hiding in the Julian Century, why the naive port was 50x slower, and what actually closed the gap.

NREL publishes a Solar Position Algorithm that computes the sun’s position for any point on Earth, any date from year -2000 to 6000, accurate to ±0.0003°. The reference implementation is C. I wanted it in JavaScript.

I assumed this was a transcription job. It took weeks, and almost all of that was one class of bug.

Both languages use IEEE 754 doubles, and that isn’t enough

JavaScript’s Number is a 64-bit IEEE 754 double. So is C’s double. People conclude from that the arithmetic is identical, and it usually is, right up until it isn’t.

C compilers are permitted to keep intermediates in wider registers, and historically x87 did exactly that at 80-bit extended precision. Optimization flags can reassociate floating-point expressions, which is invalid under IEEE semantics but common. So the “reference” output you’re matching against may itself carry extra precision that the C standard never promised you.

That means you can’t validate by diffing intermediates. You have to find the places where the algorithm is structurally precision-hostile and handle them deliberately.

The Julian Century is where the bits die

Here’s the function that cost me the most time. It’s four tokens long:

function julian_century(jd) {
    return (jd - 2451545.0) / 36525.0;
}

jd for any modern date is around 2,460,000. You’re subtracting 2,451,545 from it. Two large, nearly equal numbers.

A double carries about 15-16 significant decimal digits. At magnitude 2.46e6, roughly seven of those digits are spent on the integer part before you reach the fractional seconds you care about. The subtraction is exact in the sense that it introduces no new error, but it drops the magnitude to ~1e4 while preserving the absolute error already present. Your relative precision just fell off a cliff.

This is catastrophic cancellation, and it’s textbook. What makes it dangerous here is what happens next: jce feeds every downstream term. Nutation, obliquity, aberration, and the whole stack of periodic terms are polynomials in jce. An error introduced at this line doesn’t stay put, it propagates through several hundred subsequent operations and compounds.

The fix is to stop building a large number and then subtracting most of it away. Carry the day count and the fractional day separately for as long as possible, so the fraction keeps full mantissa precision instead of living in the low bits of a seven-digit integer. Same math, different associativity, and the associativity is the whole thing.

The other trap: truncation isn’t truncation

julian_day = integer(365.25 * (year + 4716.0))
           + integer(30.6001 * (month + 1))
           + day_decimal - 1524.5;

Note integer() rather than Math.floor or a bare cast. C’s (int) cast truncates toward zero. Math.floor rounds toward negative infinity. For positive values they agree, which is why this bug survives every test you write against modern dates and then fails somewhere in the BC range the algorithm claims to support.

The magic constants matter too. 30.6001 is not an approximation of anything physical, it’s chosen so that integer truncation lands on the right month offset. Change the truncation semantics and the constant no longer does its job.

Performance: 50x, then 3x

The naive port ran about 50x slower than C. The algorithm is a long chain of trigonometric evaluations over hundreds of periodic terms, which is close to the worst case for a JIT: heavy Math.sin/Math.cos traffic and lots of intermediate allocation if you structure the term tables as arrays of objects.

What closed most of the gap was layout and shape stability. Keep the periodic-term tables as flat numeric arrays rather than arrays of objects, so the engine sees uniform element kinds and skips the boxing. Keep every hot function monomorphic, always called with the same argument shapes, so call sites stay in the fast path instead of going megamorphic. Precompute wherever the algorithm’s structure allows a value to be hoisted out of the term loop.

That got it within 3x of C. A WebAssembly build closed it to about 1.5x, which is roughly the floor for this shape of numeric work once you’re paying the JS/WASM boundary cost at all.

The 3x number is the one I’d emphasize, though. Most numeric JavaScript that’s “too slow” isn’t slow because JavaScript is slow. It’s slow because the data layout forced the engine out of its fast paths, and that’s fixable without leaving the language.

Validation

Matched against NREL’s published test vectors and independently against ephemeris data, holding ±0.0003° on zenith, azimuth, and incidence. That precision is what makes it usable for solar positioning, daylighting analysis, and astronomy rather than just approximately right.

It’s on npm with no runtime dependencies, and it’s one of those small libraries a specific group of people genuinely needs. Those are my favorite things to ship.

The layout lesson generalized further than I expected. Performance in a mature runtime is usually determined by memory access patterns rather than by algorithmic complexity, and that same principle is what makes PHP 7’s packed hashtables a 2x win over the pointer-chasing implementation they replaced. Two different languages, a decade apart, same underlying constraint: the cache line doesn’t care what you meant.

It’s also why I reach for flat numeric arrays by default in anything with a hot loop, including the ranked lists Cascade fuses during retrieval. The instinct came from this project.