37 lines
2.2 KiB
TypeScript
37 lines
2.2 KiB
TypeScript
// Independent NOAA/Meeus solar calculator, used only by verification, never by the app.
|
|
// Equations: https://gml.noaa.gov/grad/solcalc/calcdetails.html
|
|
// https://gml.noaa.gov/grad/solcalc/main.js
|
|
// Geocentric solar coordinates + NOAA piecewise atmospheric refraction.
|
|
export function noaaPosition(date: Date, latitude: number, longitude: number) {
|
|
const rad = Math.PI / 180;
|
|
const sin = (d: number) => Math.sin(d * rad);
|
|
const cos = (d: number) => Math.cos(d * rad);
|
|
const tan = (d: number) => Math.tan(d * rad);
|
|
const t = (date.getTime() / 86400000 + 2440587.5 - 2451545) / 36525;
|
|
const l = ((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 = sin(m) * (1.914602 - t * (0.004817 + 0.000014 * t))
|
|
+ sin(2 * m) * (0.019993 - 0.000101 * t) + sin(3 * m) * 0.000289;
|
|
const omega = 125.04 - 1934.136 * t;
|
|
const lambda = l + c - 0.00569 - 0.00478 * sin(omega);
|
|
const epsilon = 23 + (26 + (21.448 - t * (46.815 + t * (0.00059 - t * 0.001813))) / 60) / 60 + 0.00256 * cos(omega);
|
|
const decl = Math.asin(sin(epsilon) * sin(lambda)) / rad;
|
|
const y = tan(epsilon / 2) ** 2;
|
|
const eq = 4 / rad * (y * sin(2 * l) - 2 * e * sin(m) + 4 * e * y * sin(m) * cos(2 * l)
|
|
- 0.5 * y ** 2 * sin(4 * l) - 1.25 * e ** 2 * sin(2 * m));
|
|
const utcMinutes = date.getUTCHours() * 60 + date.getUTCMinutes() + date.getUTCSeconds() / 60;
|
|
const solarMinutes = ((utcMinutes + eq + 4 * longitude) % 1440 + 1440) % 1440;
|
|
const ha = solarMinutes / 4 - 180;
|
|
const geometric = Math.asin(sin(latitude) * sin(decl) + cos(latitude) * cos(decl) * cos(ha)) / rad;
|
|
const azimuth = (Math.atan2(sin(ha), cos(ha) * sin(latitude) - tan(decl) * cos(latitude)) / rad + 180) % 360;
|
|
const te = tan(geometric);
|
|
let refraction = 0;
|
|
if (geometric <= 85) {
|
|
if (geometric > 5) refraction = 58.1 / te - 0.07 / te ** 3 + 0.000086 / te ** 5;
|
|
else if (geometric > -0.575) refraction = 1735 + geometric * (-518.2 + geometric * (103.4 + geometric * (-12.79 + geometric * 0.711)));
|
|
else refraction = -20.774 / te;
|
|
}
|
|
return { azimuth, altitude: geometric + refraction / 3600 };
|
|
}
|