# Eclipse solver

Apr 13th, 2024
908
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. // -------------------------------------------------------------------------- //
2. // BUILT IN MATH FUNCTIONS AND GENERAL TRIG FUNCTIONS ----------------------- //
3. // -------------------------------------------------------------------------- //
4.
5. const { PI, abs, sqrt, sign, sin, cos, tan, asin, acos, atan } = Math;
6.
7. const TAU = PI * 2;
8. const D90 = PI / 2;
9. const D30 = PI / 6;
10.
11. const toRad = (deg) => deg / 180 * PI;
13.
14. const oneWayAngleDiff = (a, b) => {
15.     return a < b ? a + TAU - b : a - b;
16. };
17.
18. // -------------------------------------------------------------------------- //
19. // CONSTANTS ---------------------------------------------------------------- //
20. // -------------------------------------------------------------------------- //
21.
22. const SEC = 1000;
23. const MIN = SEC * 60;
24. const HOUR = MIN * 60;
25. const KM = 1;
26. const MILE = 1.609344 * KM;
27. const AU = 150e6 * KM;
28. const KMPH = KM / HOUR;
29. const MPH = MILE / HOUR;
30.
31. const moonRadius = 1737.4 * KM;
32. const earthRadius = 6371 * KM;
33. const sunRadius = (1392700 / 2) * KM;
34.
35. const overrideMoonDist = 359829 * KM;
36. const overrideSunDist = 149820729 * KM;
37.
38. // -------------------------------------------------------------------------- //
39. // AUXILIAR FUNCTIONS ------------------------------------------------------- //
40. // -------------------------------------------------------------------------- //
41.
42. const parseAngle = (string) => {
43.     const sepRegex = /^\s*[°'"]\s*|\s+/;
44.     const numRegex = /^\d+(\.\d+)?/;
45.     const unitMap = {
46.         '°': 1,
47.         "'": 1 / 60,
48.         '"': 1 / 60 / 60,
49.     };
50.     string = string.trim();
51.     let sign = 1;
52.     if (string.startsWith('-')) {
53.         sign = -1;
54.         string = string.substring(1).trim();
55.     }
56.     let sum = 0;
57.     let unit = 1;
58.     while (string.length > 0) {
59.         const num = string.match(numRegex)?.[0];
60.         if (!num) {
61.             return NaN;
62.         }
63.         string = string.substring(num.length);
64.         const sep = string.match(sepRegex)?.[0];
65.         if (sep) {
66.             string = string.substring(sep.length);
67.             const short = sep.trim();
68.             if (short !== '') {
69.                 unit = unitMap[short];
70.             }
71.         } else if (string !== '') {
72.             return NaN;
73.         }
74.         sum += num * unit;
75.         unit /= 60;
76.     }
78. };
79.
80. const round = (value, n = 0) => {
81.     return Number(value.toFixed(n));
82. };
83.
84. const strAngle = (angle, pos = '', neg = '-') => {
85.     const tSec = round(toDeg(abs(angle)) * 3600, 1);
86.     const sec = round(tSec % 60, 1);
87.     const tMin = round((tSec - sec)/60);
88.     const min = round(tMin % 60);
89.     const deg = round((tMin - min)/60);
90.     const sign = angle < 0 ? neg : pos;
91.     return `\${sign}\${deg}°\${min}'\${sec}" / \${round(toDeg(angle), 6).toString()}°`;
92. };
93.
94. const numberWithNDigits = (num, n, suffix = '') => {
95.     if (abs(num) >= (10 ** n)) {
96.         return Math.round(num).toString() + suffix;
97.     }
98.     return Number(num.toPrecision(n)) + suffix;
99. };
100.
101. const strDist = (dist, n = 5) => {
102.     const km = numberWithNDigits(dist/KM, n, ' km');
103.     const mi = numberWithNDigits(dist/MILE, n, ' mi');
104.     return  `\${km} / \${mi}`;
105. };
106.
107. const strSpeed = (speed, n = 5) => {
108.     const kmph = numberWithNDigits(speed / KMPH, n, ' km/h');
109.     const mph = numberWithNDigits(speed / MPH, n, ' mph');
110.     return `\${kmph} / \${mph}`;
111. };
112.
113. const strTime = (timestamp) => {
114.     const date = new Date(timestamp).toISOString();
115.     return date.replace(/^.*T/, '').replace('Z', ' (UTC)');
116. };
117.
118. const strDuration = (duration) => {
119.     const tSec = round(duration / SEC);
120.     const sec = tSec % 60;
121.     const tMin = (tSec - sec)/60;
122.     const min = tMin % 60;
123.     const hrs = (tMin - min)/60;
124.     const strSec = sec.toString().padStart(2, '0');
125.     const strMin = min.toString().padStart(2, '0');
126.     const strHrs = hrs.toString().padStart(2, '0');
127.     return `\${strHrs}h\${strMin}m\${strSec}s`;
128. };
129.
130. // -------------------------------------------------------------------------- //
131. // DATA PARSING ------------------------------------------------------------- //
132. // -------------------------------------------------------------------------- //
133.
134. const table = [];
135.
136. const addAlmanacData = (text) => {
137.     let date = null;
138.     const lines = text.trim().split(/\n/);
139.     for (let line of lines) {
140.         line = line.trim();
141.         if (line.startsWith('#')) {
142.             continue;
143.         }
144.         if (/^date:/i.test(line)) {
145.             const strDate = line.replace(/^date:/i, '').trim();
146.             const isoDate = strDate + 'T00:00:00Z';
147.             date = new Date(isoDate).getTime();
148.             continue;
149.         }
150.         const [ hour, ...angles ] = line.trim().split(/\s*\|\s*/);
151.         const time = date + hour * HOUR;
152.         const [ sunGHA, sunDec, moonGHA, moonDec, moonHP ] = angles.map(parseAngle);
153.         const obj = { time, sunGHA, sunDec, moonGHA, moonDec, moonHP };
154.         table.push(obj);
155.     }
156. };
157.
158. const interpolate = (time, field) => {
159.     for (let i=1; i<table.length; i++) {
160.         const bef = table[i - 1];
161.         if (time < bef.time) {
162.             continue;
163.         }
164.         const aft = table[i];
165.         if (time > aft.time) {
166.             continue;
167.         }
168.         const t = (time - bef.time) / (aft.time - bef.time);
169.         const a = bef[field];
170.         const b = aft[field];
171.         return a + (b - a) * t;
172.     }
173.     throw new Error(`time \${new Date(time).toISOString()} is out of bounds`);
174. };
175.
176. // -------------------------------------------------------------------------- //
177. // VECTOR OPERATIONS (LINEAR ALGEBRA) --------------------------------------- //
178. // -------------------------------------------------------------------------- //
179.
180. const subVec = ([ ax, ay, az ], [ bx, by, bz ]) => {
181.     return [ ax - bx, ay - by, az - bz ];
182. };
183.
184. const sumVec = ([ ax, ay, az ], [ bx, by, bz ]) => {
185.     return [ ax + bx, ay + by, az + bz ];
186. };
187.
188. const normalizeVec = ([ x, y, z ]) => {
189.     const len = sqrt(x*x + y*y + z*z);
190.     return [ x/len, y/len, z/len ];
191. };
192.
193. const dotProd = (a, b) => {
194.     const [ ax, ay, az ] = a;
195.     const [ bx, by, bz ] = b;
196.     return (ax * bx) + (ay * by) + (az * bz);
197. };
198.
199. const scaleVec = (vec, val) => {
200.     const [ x, y, z ] = vec;
201.     return [ x*val, y*val, z*val ];
202. };
203.
204. const vecLen = ([ x, y, z ]) => {
205.     return sqrt(x*x + y*y + z*z);
206. };
207.
208. const rotVecY = ([ x, y, z ], angle) => {
209.     const s = sin(angle);
210.     const c = cos(angle);
211.     return [ x*c - z*s, y, z*c + x*s ];
212. };
213.
214. const rotVecX = ([ x, y, z ], angle) => {
215.     const s = sin(angle);
216.     const c = cos(angle);
217.     return [ x, y*c + z*s, z*c - y*s ];
218. };
219.
220. const rotVecZ = ([ x, y, z ], angle) => {
221.     const s = sin(angle);
222.     const c = cos(angle);
223.     return [ x*c + y*s, y*c - x*s, z ];
224. };
225.
226. const vecIsValid = ([ x, y, z ]) => {
227.     return !(isNaN(x) || isNaN(y) || isNaN(z));
228. };
229.
230. // -------------------------------------------------------------------------- //
231. // LAT/LONG MATH ------------------------------------------------------------ //
232. // -------------------------------------------------------------------------- //
233.
234. // Geocentric
235. const GC = {
236.     vecToGP: (vec) => {
237.         const [ x, y, z ] = normalizeVec(vec);
238.         const lat = asin(z);
239.         const len = sqrt(x*x + y*y);
240.         const lon = len === 0 ? 0 : acos(x/len) * (y < 0 ? -1 : 1);
241.         return [ lat, lon ];
242.     },
243.     gpToUnitVec: (gp) => {
244.         const [ lat, lon ] = gp;
245.         const x = + cos(lat) * cos(lon);
246.         const y = + cos(lat) * sin(lon);
247.         const z = + sin(lat);
248.         return [ x, y, z ];
249.     },
250.     gpDistToVec: (gp, dist) => {
251.         const vec = GC.gpToUnitVec(gp);
252.         return scaleVec(vec, dist);
253.     },
254.     calcAngle: (gp1, gp2) => {
255.         const aVec = GC.gpToUnitVec(gp1);
256.         const bVec = GC.gpToUnitVec(gp2);
257.         return acos(dotProd(aVec, bVec));
258.     },
259.     calcDistance: (gp1, gp2) => {
260.         return GC.calcAngle(gp1, gp2) * earthRadius;
261.     },
262.     calcAzimuth: (gp1, gp2) => {
263.         const [ lat1, lon1 ] = gp1;
264.         const [ lat2, lon2 ] = gp2;
265.         const dlon = lon2 - lon1;
266.         const x = cos(lat2) * cos(dlon);
267.         const y = cos(lat2) * sin(dlon);
268.         const z = sin(lat2);
269.         const newZ = z*cos(lat1) - x*sin(lat1);
270.         const east = acos(newZ/sqrt(newZ ** 2 + y ** 2));
271.         return y >= 0 ? east : TAU - east;
272.     },
273. };
274.
275. // -------------------------------------------------------------------------- //
276. // ECLIPSE MATH ------------------------------------------------------------- //
277. // -------------------------------------------------------------------------- //
278.
279. const ghaToLon = (gha) => {
280.     return (TAU*1.5 - gha) % TAU - PI;
281. };
282.
283. const getMoonGhaDec = (time) => {
284.     const gha = interpolate(time, 'moonGHA');
285.     const dec = interpolate(time, 'moonDec');
286.     return [ gha, dec ];
287. };
288.
289. const getSunGhaDec = (time) => {
290.     const gha = interpolate(time, 'sunGHA');
291.     const dec = interpolate(time, 'sunDec');
292.     return [ gha, dec ];
293. };
294.
295. const hpToDist = (hp) => {
296.     const dist = earthRadius / tan(hp);
297.     return dist;
298. };
299.
300. const findGHAMeridianMatch = () => {
301.     for (let i=1; i<table.length; ++i) {
302.         const bef = table[i - 1];
303.         const aft = table[i];
304.         const sunApproaching = oneWayAngleDiff(bef.moonGHA, bef.sunGHA) < D30;
305.         const sunPasses = oneWayAngleDiff(aft.sunGHA, aft.moonGHA) < D30;
306.         if (!sunApproaching || !sunPasses) {
307.             continue;
308.         }
309.         const dt = aft.time - bef.time;
310.         const dSun = aft.sunGHA - bef.sunGHA;
311.         const dMoon = aft.moonGHA - bef.moonGHA;
312.         const t = (bef.moonGHA - bef.sunGHA) / (dSun - dMoon);
313.         const time = bef.time + dt * t;
314.         return time;
315.     }
316.     return null;
317. };
318.
319. const ghaDecDistToVec = (gha, dec, dist) => {
320.     const gp = [ dec, ghaToLon(gha) ];
321.     return GC.gpDistToVec(gp, dist);
322. };
323.
324. const getMoonDist = (time) => {
325.     return overrideMoonDist ?? hpToDist(interpolate(time, 'moonHP'));
326. };
327.
328. const getMoonVec = (time) => {
329.     const [ gha, dec ] = getMoonGhaDec(time);
330.     const dist = getMoonDist(time);
332. };
333.
334. const getSunDist = (time) => {
335.     return overrideSunDist ?? AU;
336. };
337.
338. const getSunVec = (time) => {
339.     const [ gha, dec ] = getSunGhaDec(time);
340.     const dist = getSunDist(time);
342. };
343.
344. const getCenterShadowRay = (moonVec, sunVec) => {
345.     const dir = normalizeVec(subVec(moonVec, sunVec));
346.     const origin = moonVec;
347.     return { origin, dir };
348. };
349.
350. const rayEarthIntersection = ({ origin, dir }) => {
351.     const midT = dotProd(subVec([ 0, 0, 0 ], origin), dir);
352.     if (midT < 0) {
353.         return [ NaN, NaN, NaN ];
354.     }
355.     const midVec = sumVec(origin, scaleVec(dir, midT));
356.     const centerDist = vecLen(midVec);
357.     if (centerDist > earthRadius) {
358.         return [ NaN, NaN, NaN ];
359.     }
360.     const tDiff = sqrt(earthRadius**2 - centerDist**2);
361.     const t = midT - tDiff;
362.     return sumVec(origin, scaleVec(dir, t));
363. };
364.
365. const calcUmbraAngle = (moonSunDist) => {
367. };
368.
369. const calcPenumbraAngle = (moonSunDist) => {
371. };
372.
373. const buildShadowEdgeRay = (moonVec, sunVec, azimuth, calcAngle) => {
374.     const moonSunDist = vecLen(subVec(moonVec, sunVec));
375.     const angle = calcAngle(moonSunDist);
376.
377.     let origin = [ 0, 0, moonRadius ];
378.     let dir = [ -1, 0, 0 ];
379.
380.     origin = rotVecY(origin, angle);
381.     dir = rotVecY(dir, angle)
382.
383.     origin = rotVecX(origin, azimuth);
384.     dir = rotVecX(dir, azimuth)
385.
386.     const [ lat, lon ] = GC.vecToGP(subVec(sunVec, moonVec));
387.
388.     origin = rotVecY(origin, lat);
389.     dir = rotVecY(dir, lat);
390.
391.     origin = rotVecZ(origin, -lon);
392.     dir = rotVecZ(dir, -lon);
393.
394.     origin = sumVec(origin, moonVec);
395.
396.     return { origin, dir };
397. };
398.
399. const getShadowEdgeCoord = (moonVec, sunVec, azimuth, calcAngle) => {
400.     return GC.vecToGP(
401.         rayEarthIntersection(
403.         ),
404.     );
405. };
406.
407. const calcShadowSize = (moonVec, sunVec, azimuth, calcAngle) => {
408.     const az1 = azimuth;
409.     const az2 = (azimuth + PI) % TAU;
410.     const a = getShadowEdgeCoord(moonVec, sunVec, az1, calcAngle);
411.     const b = getShadowEdgeCoord(moonVec, sunVec, az2, calcAngle);
412.     return GC.calcDistance(a, b);
413. };
414.
415. const calcVecs = (time) => {
416.     const moonVec = getMoonVec(time);
417.     const [ sunGHA, sunDec ] = getSunGhaDec(time);
418.     const sunGP = [ sunDec, ghaToLon(sunGHA) ];
419.     const sunDist = getSunDist(time);
420.     const sunVec = ghaDecDistToVec(sunGHA, sunDec, sunDist);
421.     const locVec = rayEarthIntersection(getCenterShadowRay(moonVec, sunVec));
422.     const loc = GC.vecToGP(locVec);
423.     return { sunGP, moonVec, sunVec, loc, locVec };
424. };
425.
426. const hasTotality = (locVec, time) => {
427.     const { sunVec, moonVec } = calcVecs(time);
430.     const sunDir = normalizeVec(subVec(sunVec, locVec));
431.     const moonDir = normalizeVec(subVec(moonVec, locVec));
432.     const prod = dotProd(moonDir, sunDir);
433.     if (prod >= 1) {
434.         return true;
435.     }
436.     const angDist = acos(dotProd(moonDir, sunDir));
438.     return angDist <= gap;
439. };
440.
442.     const dist = vecLen(subVec(locVec, bodyVec));
444. };
445.
446. const findTotalityEdgeAt = (locVec, startTime, endTime) => {
447.     let at = startTime;
448.     let av = hasTotality(locVec, at);
449.     let bt = endTime;
450.     let bv = hasTotality(locVec, bt);
451.     if (av === bv) {
452.         return null;
453.     }
454.     while (bt - at > 1) {
455.         const mt = (at + bt)/2;
456.         const mv = hasTotality(locVec, mt);
457.         if (mv === av) {
458.             at = mt;
459.             av = mv;
460.         } else {
461.             bt = mt;
462.             bv = mv;
463.         }
464.     }
465.     return (at + bt)/2;
466. };
467.
468. // -------------------------------------------------------------------------- //
469. // MAIN --------------------------------------------------------------------- //
470. // -------------------------------------------------------------------------- //
471.
473.     Date: 2024-04-08
474.     # Hour | Sun GHA    | Sun Dec  | Moon GHA   | Moon Dec | Moon HP
475.     #------|------------|----------|------------|----------|---------
476.       18   |  89° 35.5' | 7° 35.2' |  89° 54.4' | 7° 48.9' | 60.9'
477.       19   | 104° 35.6' | 7° 36.2' | 104° 23.2' | 8° 06.3' | 60.9'
478. `);
479.
480. const peakTime = findGHAMeridianMatch();
481. const { sunGP, moonVec, sunVec, loc, locVec } = calcVecs(peakTime);
482. const [ lat, lon ] = loc;
483. const umbraSizeNS = calcShadowSize(moonVec, sunVec,   0, calcUmbraAngle);
484. const umbraSizeEW = calcShadowSize(moonVec, sunVec, D90, calcUmbraAngle);
485. const penumbraSizeNS = calcShadowSize(moonVec, sunVec,   0, calcPenumbraAngle);
486. const penumbraSizeEW = calcShadowSize(moonVec, sunVec, D90, calcPenumbraAngle);
487. const sunAlt = D90 - GC.calcAngle(sunGP, loc);
488.
489. const dtTime = 1 * SEC;
490. const nextTime = peakTime + dtTime;
491. const nextLoc = calcVecs(nextTime).loc;
492. const dtDist = GC.calcDistance(loc, nextLoc);
493. const speed = dtDist/dtTime;
494. const dir = GC.calcAzimuth(loc, nextLoc);
495. const pathWidth = calcShadowSize(moonVec, sunVec, dir + D90, calcUmbraAngle);
496. const startTime = findTotalityEdgeAt(locVec, table[0].time, peakTime);
497. const endTime = findTotalityEdgeAt(locVec, peakTime, table.at(-1).time);
498. const duration = endTime - startTime;
501.
502. console.log(' 1. Latitude:', strAngle(lat, 'N ', 'S '));
503. console.log(' 2. Longitude:', strAngle(lon, 'E ', 'W '));
504. console.log(' 3. Time:', strTime(peakTime));
505. console.log(' 4. Sun Alt.:', strAngle(sunAlt));