Resources Numerical accuracy Calculation policy

Sun position

Works out where the sun is for a date, time and place: azimuth, altitude with refraction, declination and the equation of time.

How the precision is kept

64-bit arithmetic on published series. The accuracy is that of the series and its constants, which are written to the digits the source publishes, not of the arithmetic. The file comment gives about a minute of arc for the years 1800 to 2100.

Calculations under this policy

3 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.

drive.js CivilCraft 0.8.1 64-bit 3 calculations Reviewed Nov 5, 2026

Tests bin/test-civilcraft-drive.mjs, bin/test-civilcraft-massing.mjs

Sun azimuth and altitude

64-bit

The sun's direction and height above the horizon at a moment and a place.

JD = ms ÷ 86 400 000 + 2 440 587.5; T = (JD − 2 451 545) ÷ 36 525; cos Z = sin φ·sin δ + cos φ·cos δ·cos H

The Julian day is a number near 2.46 million, so a 64-bit value resolves it to about 40 microseconds; the century fraction T is then a difference of two such numbers. The series constants are written as published (for example 280.46646 and 36000.76983). Both acos arguments are clamped to [−1, 1].

drive.js
lines 50 to 77
CivilCraft 0.8.1

export function sunPosition(date, lat, lon) {
  const ms = date instanceof Date ? date.getTime() : Number(date);
  const jd = ms / 86400000 + 2440587.5, T = (jd - 2451545) / 36525;
  const L0 = (((280.46646 + T * (36000.76983 + T * 0.0003032)) % 360) + 360) % 360;
  const M = 357.52911 + T * (35999.05029 - 0.0001537 * T);
  const e = 0.016708634 - T * (0.000042037 + 0.0000001267 * T);
  const C = Math.sin(M * R) * (1.914602 - T * (0.004817 + 0.000014 * T)) + Math.sin(2 * M * R) * (0.019993 - 0.000101 * T) + Math.sin(3 * M * R) * 0.000289;
  const omega = 125.04 - 1934.136 * T, lambda = L0 + C - 0.00569 - 0.00478 * Math.sin(omega * R);
  const eps0 = 23 + (26 + (21.448 - T * (46.815 + T * (0.00059 - T * 0.001813))) / 60) / 60, eps = eps0 + 0.00256 * Math.cos(omega * R);
  const decl = Math.asin(Math.sin(eps * R) * Math.sin(lambda * R)) / R;
  const y = Math.tan((eps / 2) * R) ** 2;
  const eot = (4 / R) * (y * Math.sin(2 * L0 * R) - 2 * e * Math.sin(M * R) + 4 * e * y * Math.sin(M * R) * Math.cos(2 * L0 * R) - 0.5 * y * y * Math.sin(4 * L0 * R) - 1.25 * e * e * Math.sin(2 * M * R));
  const d = new Date(ms), minutes = d.getUTCHours() * 60 + d.getUTCMinutes() + d.getUTCSeconds() / 60 + d.getUTCMilliseconds() / 60000;
  let tst = (minutes + eot + 4 * lon) % 1440; if (tst < 0) tst += 1440;
  const ha = tst / 4 < 0 ? tst / 4 + 180 : tst / 4 - 180;
  const cosZ = Math.sin(lat * R) * Math.sin(decl * R) + Math.cos(lat * R) * Math.cos(decl * R) * Math.cos(ha * R);
  const zen = Math.acos(Math.max(-1, Math.min(1, cosZ))) / R, elev = 90 - zen;
  // azimuth from north, clockwise
  const sz = Math.sin(zen * R);
  let az;
  if (Math.abs(sz) < 1e-9) az = lat > decl ? 180 : 0;
  else {
    const c = (Math.sin(lat * R) * Math.cos(zen * R) - Math.sin(decl * R)) / (Math.cos(lat * R) * sz);
    const a = Math.acos(Math.max(-1, Math.min(1, c))) / R;
    az = ha > 0 ? (a + 180) % 360 : (540 - a) % 360;
  }
  return { azimuth: az, altitude: elev + refraction(elev), declination: decl, eot };
}

Atmospheric refraction

64-bit

The lift of the sun's apparent altitude by the atmosphere, added to the geometric altitude.

h > 5°: 58.1 ÷ tan h − 0.07 ÷ tan³h + 0.000086 ÷ tan⁵h seconds of arc

Three empirical formulas chosen by altitude, with constants of two to four digits. They do not join exactly at 5° and at −0.575°, so the correction steps slightly there.

drive.js
lines 78 to 87
CivilCraft 0.8.1

/** NOAA's approximation of atmospheric refraction (degrees) at an apparent elevation. */
function refraction(h) {
  if (h > 85) return 0;
  const te = Math.tan(h * R);
  let r;
  if (h > 5) r = 58.1 / te - 0.07 / te ** 3 + 0.000086 / te ** 5;
  else if (h > -0.575) r = 1735 + h * (-518.2 + h * (103.4 + h * (-12.79 + h * 0.711)));
  else r = -20.774 / te;
  return r / 3600;
}

Unit vector toward the sun

64-bit

The sun's direction as east, north and up components.

(cos h·sin a, cos h·cos a, sin h)

Three products of sines and cosines of the 64-bit angles. The vector has unit length to rounding and is not renormalized.

drive.js
lines 88 to 91
CivilCraft 0.8.1

export function sunVector({ azimuth, altitude }) {
  const a = azimuth * R, h = altitude * R;
  return [Math.cos(h) * Math.sin(a), Math.cos(h) * Math.cos(a), Math.sin(h)];
}

In detail

Method: The file comment names NOAA's solar calculator equations (after Meeus): mean longitude and anomaly, equation of center, apparent longitude, obliquity, declination, equation of time, hour angle, then zenith and azimuth; NOAA's refraction approximation in three ranges of altitude.
Number types: 64-bit numbers. Angles in degrees, converted with π ÷ 180 at each trigonometric call.
Tolerances: Arguments of acos are clamped to [−1, 1]. With the sun within about 1e-9 of the zenith (sine of the zenith angle below 1e-9) the azimuth is set to 180° or 0°. No refraction above 85° of altitude.
Limits: A truncated series: about a minute of arc for 1800 to 2100 according to the file comment, which was not re-derived in this review. The refraction formula is for a standard atmosphere.

Approach: 64-bit, stated limit

The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.

Every policy that takes this approach