Calcuate direction to sun based on time of day, latitude, axial tilt

sky.js

const { cos, PI, sin } = Math;

let time = 5; // 0-24
const LAT = 45 * PI / 180;
const AXIAL_TILT = 23.5 * PI / 180;
const SOLAR_ORBIT = PI/2; // 0 - tilted towards the sun = summer solstice, PI = winter solstice

export function update() {

  //////////////////////////////////////////////////////////////////////////
  // Calculate angle to sun
  let cosv, sinv, sn0, sn1, sn2; // eslint-disable-line one-var
  let sun_dir = vec3(0, 0, -1); // midnight, with no tilt, at the equator, sun is straight down
  // rotate around the sun (rotate counterclockwise? around Y)
  cosv = cos(SOLAR_ORBIT);
  sinv = sin(SOLAR_ORBIT);
  sn0 = sun_dir[0];
  sn2 = sun_dir[2];
  sun_dir[0] = cosv * sn0 + sinv * sn2;
  sun_dir[2] = -sinv * sn0 + cosv * sn2;
  // apply axial tilt (rotate counterclockwise? around X)
  cosv = cos(AXIAL_TILT);
  sinv = sin(AXIAL_TILT);
  sn1 = sun_dir[1];
  sn2 = sun_dir[2];
  sun_dir[1] = cosv * sn1 - sinv * sn2;
  sun_dir[2] = sinv * sn1 + cosv * sn2;
  // rotate clockwise around Y by time
  let time_angle = time / 24 * PI * 2 + SOLAR_ORBIT;
  cosv = cos(-time_angle);
  sinv = sin(-time_angle);
  sn0 = sun_dir[0];
  sn2 = sun_dir[2];
  sun_dir[0] = cosv * sn0 + sinv * sn2;
  sun_dir[2] = -sinv * sn0 + cosv * sn2;
  // rotate counterclockwise by latitude around X
  cosv = cos(LAT);
  sinv = sin(LAT);
  sn1 = sun_dir[1];
  sn2 = sun_dir[2];
  sun_dir[1] = cosv * sn1 - sinv * sn2;
  sun_dir[2] = sinv * sn1 + cosv * sn2;
  // We now have a directional vector pointing towards the sun
  v3copy(sky_props.sun_pos, sun_dir);
  v3scale(sky_props.light_dir_ws, sun_dir, -1);
}