summaryrefslogtreecommitdiff
path: root/tools/sky/astro_test.js
blob: b299d163190e48e037bc34835aef98ecef49740b (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
// Checks public/js/astro.js against JPL Horizons and Meeus.
//
//   deno test tools/sky/astro_test.js
import '../../public/js/astro.js';

const { Astro } = globalThis;
const RAD = Math.PI / 180;
const gap = (a, b) => Math.abs((((a - b + 180) % 360) + 360) % 360 - 180);
const ok = (cond, msg) => { if (!cond) throw new Error(msg); };

function eclipticOf(dir, ms) {
  // Inverse of Astro.equatorial, to compare with Horizons' ecliptic longitudes.
  const ob = (23.4393 - 3.563e-7 * (Astro.julianDay(ms) - 2451543.5)) * RAD;
  const [x, y, z] = dir;
  const ye = y * Math.cos(ob) + z * Math.sin(ob);
  const ze = -y * Math.sin(ob) + z * Math.cos(ob);
  return { lon: ((Math.atan2(ye, x) / RAD) + 360) % 360, lat: Math.asin(ze) / RAD };
}

Deno.test('planets are where Horizons puts them (geocentric, of date)', () => {
  // From android-toolkit's AstroTest: apparent ecliptic longitude and latitude of the date.
  const refs = {
    '2026-10-31T20:00:00Z': { sun: [218.4102216, 0.0001252], moon: [115.1119618, 2.6704709], mercury: [226.8441830, -1.5976931], venus: [206.3301706, -4.9067015], mars: [138.5376903, 1.5969459], jupiter: [144.3028948, 0.6895096], saturn: [9.2969189, -2.6851943] },
    '2028-08-17T06:00:00Z': { sun: [144.8636375, 0.0000135], moon: [101.0244024, 1.0776734], mercury: [165.1873864, 0.4638201], venus: [99.2952316, -2.9103916], mars: [108.2009535, 0.7458269], jupiter: [178.6149976, 1.1269943], saturn: [41.2847987, -2.5175988] },
  };
  for (const [iso, bodies] of Object.entries(refs)) {
    const ms = Date.parse(iso);
    for (const [name, [lon, lat]] of Object.entries(bodies)) {
      const e = Astro.ecliptic(name, ms);
      ok(gap(e.lon, lon) < 0.05, `${name} ${iso} lon ${e.lon} vs ${lon}`);
      ok(Math.abs(e.lat - lat) < 0.05, `${name} ${iso} lat ${e.lat} vs ${lat}`);
      const back = eclipticOf(Astro.equatorial(e.lon, e.lat, ms), ms);
      ok(gap(back.lon, e.lon) < 1e-6 && Math.abs(back.lat - e.lat) < 1e-6, `${name} round trip`);
    }
  }
});

Deno.test("Seattle's sky matches Horizons (airless azimuth, elevation, magnitude)", () => {
  // Horizons, observer at 47.6062 N, 122.3321 W, 2026-11-01 00:00 UT (JD 2461345.5).
  const ms = Date.parse('2026-11-01T00:00:00Z');
  const refs = {
    sun: [239.356839801, 7.223273431],
    moon: [338.838357584, -17.481568478, -10.609],
    mercury: [231.106696194, 8.172959876, 3.268],
    venus: [249.495509600, -1.327997551, -4.170],
    mars: [316.936325646, -14.742817319, 0.789],
    jupiter: [310.975220809, -14.507738850, -2.003],
    saturn: [87.986528867, -0.198109819, 0.462],
    uranus: [33.812741373, -14.521612899, 5.611],
  };
  const horizon = Astro.horizonMatrix(ms, 47.6062, -122.3321);
  const sky = Astro.bodies(ms, 47.6062, -122.3321);
  const found = { sun: { dir: sky.sun }, moon: sky.moon };
  for (const p of sky.planets) found[p.name] = p;
  for (const [name, [az, alt, mag]] of Object.entries(refs)) {
    const h = Astro.altAz(found[name].dir, horizon);
    ok(Math.abs(h.alt - alt) < 0.05, `${name} alt ${h.alt} vs ${alt}`);
    ok(gap(h.az, az) < 0.1, `${name} az ${h.az} vs ${az}`);
    // Schlyter's magnitudes are rough for Mercury near inferior conjunction, as it is here.
    if (mag !== undefined) ok(Math.abs(found[name].magnitude - mag) < (name === 'mercury' ? 0.5 : 0.25), `${name} mag ${found[name].magnitude} vs ${mag}`);
  }
});

Deno.test('precession matches Meeus example 21.b', () => {
  // Theta Persei, J2000 to 2028 Nov 13.19 TD, with its proper motion (+0.99 s, -2.58") taken back out.
  const ms = (2462088.69 - 2440587.5) * 86400000;
  const ra = (2 + 44 / 60 + 11.986 / 3600) * 15 * RAD;
  const dec = (49 + 13 / 60 + 42.48 / 3600) * RAD;
  const v = Astro.apply(Astro.precessionMatrix(ms), [Math.cos(dec) * Math.cos(ra), Math.cos(dec) * Math.sin(ra), Math.sin(dec)]);
  const ra2 = Math.atan2(v[1], v[0]) / RAD / 15;
  const dec2 = Math.asin(v[2]) / RAD;
  const wantRa = 2 + 46 / 60 + (11.331 - 0.99) / 3600;
  const wantDec = 49 + 20 / 60 + (54.54 + 2.58) / 3600;
  ok(Math.abs(ra2 - wantRa) * 3600 * 15 < 1, `ra ${ra2} vs ${wantRa}`);
  ok(Math.abs(dec2 - wantDec) * 3600 < 1, `dec ${dec2} vs ${wantDec}`);
});