From 5602eadf68570ee99bccf69af98c23bc034eb16f Mon Sep 17 00:00:00 2001 From: Klaus-Peter Bernschneider Date: Tue, 15 Mar 2022 09:22:07 +0100 Subject: [PATCH] Add files via upload --- modules/suncalc.js | 298 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 298 insertions(+) create mode 100644 modules/suncalc.js diff --git a/modules/suncalc.js b/modules/suncalc.js new file mode 100644 index 000000000..b1af0a0d9 --- /dev/null +++ b/modules/suncalc.js @@ -0,0 +1,298 @@ +/* Module suncalc.js + (c) 2011-2015, Vladimir Agafonkin + SunCalc is a JavaScript library for calculating sun/moon position and light phases. + https://github.com/mourner/suncalc + +PB: Usage: +E.setTimeZone(2); // 1 = MEZ, 2 = MESZ +SunCalc = require("suncalc.js"); +pos = SunCalc.getPosition(Date.now(), 53.3, 10.1); +times = SunCalc.getTimes(Date.now(), 53.3, 10.1); +rise = times.sunrise; // Date object +rise_str = rise.getHours() + ':' + rise.getMinutes(); //hh:mm +*/ +var exports={}; + +// shortcuts for easier to read formulas + +var PI = Math.PI, + sin = Math.sin, + cos = Math.cos, + tan = Math.tan, + asin = Math.asin, + atan = Math.atan2, + acos = Math.acos, + rad = PI / 180; + +// sun calculations are based on http://aa.quae.nl/en/reken/zonpositie.html formulas + +// date/time constants and conversions + +var dayMs = 1000 * 60 * 60 * 24, + J1970 = 2440588, + J2000 = 2451545; + +function toJulian(date) { return date.valueOf() / dayMs - 0.5 + J1970; } +function fromJulian(j) { return new Date((j + 0.5 - J1970) * dayMs); } // PB: onece removed + 0.5; included it again 4 Jan 2021 +function toDays(date) { return toJulian(date) - J2000; } + + +// general calculations for position + +var e = rad * 23.4397; // obliquity of the Earth + +function rightAscension(l, b) { return atan(sin(l) * cos(e) - tan(b) * sin(e), cos(l)); } +function declination(l, b) { return asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l)); } + +function azimuth(H, phi, dec) { return atan(sin(H), cos(H) * sin(phi) - tan(dec) * cos(phi)); } +function altitude(H, phi, dec) { return asin(sin(phi) * sin(dec) + cos(phi) * cos(dec) * cos(H)); } + +function siderealTime(d, lw) { return rad * (280.16 + 360.9856235 * d) - lw; } + +function astroRefraction(h) { + if (h < 0) // the following formula works for positive altitudes only. + h = 0; // if h = -0.08901179 a div/0 would occur. + + // formula 16.4 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998. + // 1.02 / tan(h + 10.26 / (h + 5.10)) h in degrees, result in arc minutes -> converted to rad: + return 0.0002967 / Math.tan(h + 0.00312536 / (h + 0.08901179)); +} + +// general sun calculations + +function solarMeanAnomaly(d) { return rad * (357.5291 + 0.98560028 * d); } + +function eclipticLongitude(M) { + + var C = rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M)), // equation of center + P = rad * 102.9372; // perihelion of the Earth + + return M + C + P + PI; +} + +function sunCoords(d) { + + var M = solarMeanAnomaly(d), + L = eclipticLongitude(M); + + return { + dec: declination(L, 0), + ra: rightAscension(L, 0) + }; +} + +// calculates sun position for a given date and latitude/longitude + +exports.getPosition = function (date, lat, lng) { + + var lw = rad * -lng, + phi = rad * lat, + d = toDays(date), + + c = sunCoords(d), + H = siderealTime(d, lw) - c.ra; + + return { + azimuth: Math.round((azimuth(H, phi, c.dec) / rad + 180) % 360), // PB: converted to deg + altitude: Math.round( altitude(H, phi, c.dec) / rad) // PB: converted to deg + }; +}; + + +// sun times configuration (angle, morning name, evening name) + +var times = [ + [-0.833, 'sunrise', 'sunset' ] +]; + +// calculations for sun times +var J0 = 0.0009; + +function julianCycle(d, lw) { return Math.round(d - J0 - lw / (2 * PI)); } + +function approxTransit(Ht, lw, n) { return J0 + (Ht + lw) / (2 * PI) + n; } +function solarTransitJ(ds, M, L) { return J2000 + ds + 0.0053 * sin(M) - 0.0069 * sin(2 * L); } + +function hourAngle(h, phi, d) { return acos((sin(h) - sin(phi) * sin(d)) / (cos(phi) * cos(d))); } +function observerAngle(height) { return -2.076 * Math.sqrt(height) / 60; } + +// returns set time for the given sun altitude +function getSetJ(h, lw, phi, dec, n, M, L) { + + var w = hourAngle(h, phi, dec), + a = approxTransit(w, lw, n); + return solarTransitJ(a, M, L); +} + + +// calculates sun times for a given date, latitude/longitude, and, optionally, +// the observer height (in meters) relative to the horizon + +exports.getTimes = function (date, lat, lng, height) { + + height = height || 0; + + var lw = rad * -lng, + phi = rad * lat, + + dh = observerAngle(height), + + d = toDays(date), + n = julianCycle(d, lw), + ds = approxTransit(0, lw, n), + + M = solarMeanAnomaly(ds), + L = eclipticLongitude(M), + dec = declination(L, 0), + + Jnoon = solarTransitJ(ds, M, L), + + i, len, time, h0, Jset, Jrise; + + + var result = { + solarNoon: fromJulian(Jnoon), + nadir: fromJulian(Jnoon - 0.5) + }; + + for (i = 0, len = times.length; i < len; i += 1) { + time = times[i]; + h0 = (time[0] + dh) * rad; + + Jset = getSetJ(h0, lw, phi, dec, n, M, L); + Jrise = Jnoon - (Jset - Jnoon); + + result[time[1]] = fromJulian(Jrise); + result[time[2]] = fromJulian(Jset); + } + + return result; +}; + + +// moon calculations, based on http://aa.quae.nl/en/reken/hemelpositie.html formulas + +function moonCoords(d) { // geocentric ecliptic coordinates of the moon + + var L = rad * (218.316 + 13.176396 * d), // ecliptic longitude + M = rad * (134.963 + 13.064993 * d), // mean anomaly + F = rad * (93.272 + 13.229350 * d), // mean distance + + l = L + rad * 6.289 * sin(M), // longitude + b = rad * 5.128 * sin(F), // latitude + dt = 385001 - 20905 * cos(M); // distance to the moon in km + + return { + ra: rightAscension(l, b), + dec: declination(l, b), + dist: dt + }; +} + +getMoonPosition = function (date, lat, lng) { + + var lw = rad * -lng, + phi = rad * lat, + d = toDays(date), + + c = moonCoords(d), + H = siderealTime(d, lw) - c.ra, + h = altitude(H, phi, c.dec), + // formula 14.1 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998. + pa = atan(sin(H), tan(phi) * cos(c.dec) - sin(c.dec) * cos(H)); + + h = h + astroRefraction(h); // altitude correction for refraction + + return { + azimuth: azimuth(H, phi, c.dec), + altitude: h, + distance: c.dist, + parallacticAngle: pa + }; +}; + + +// calculations for illumination parameters of the moon, +// based on http://idlastro.gsfc.nasa.gov/ftp/pro/astro/mphase.pro formulas and +// Chapter 48 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998. + +getMoonIllumination = function (date) { + + var d = toDays(date || new Date()), + s = sunCoords(d), + m = moonCoords(d), + + sdist = 149598000, // distance from Earth to Sun in km + + phi = acos(sin(s.dec) * sin(m.dec) + cos(s.dec) * cos(m.dec) * cos(s.ra - m.ra)), + inc = atan(sdist * sin(phi), m.dist - sdist * cos(phi)), + angle = atan(cos(s.dec) * sin(s.ra - m.ra), sin(s.dec) * cos(m.dec) - + cos(s.dec) * sin(m.dec) * cos(s.ra - m.ra)); + + return { + fraction: (1 + cos(inc)) / 2, + phase: 0.5 + 0.5 * inc * (angle < 0 ? -1 : 1) / Math.PI, + angle: angle + }; +}; + + +function hoursLater(date, h) { + return new Date(date.valueOf() + h * dayMs / 24); +} + +// calculations for moon rise/set times are based on http://www.stargazing.net/kepler/moonrise.html article + +getMoonTimes = function (date, lat, lng, inUTC) { + var t = new Date(date); + if (inUTC) t.setUTCHours(0, 0, 0, 0); + else t.setHours(0, 0, 0, 0); + + var hc = 0.133 * rad, + h0 = SunCalc.getMoonPosition(t, lat, lng).altitude - hc, + h1, h2, rise, set, a, b, xe, ye, d, roots, x1, x2, dx; + + // go in 2-hour chunks, each time seeing if a 3-point quadratic curve crosses zero (which means rise or set) + for (var i = 1; i <= 24; i += 2) { + h1 = SunCalc.getMoonPosition(hoursLater(t, i), lat, lng).altitude - hc; + h2 = SunCalc.getMoonPosition(hoursLater(t, i + 1), lat, lng).altitude - hc; + + a = (h0 + h2) / 2 - h1; + b = (h2 - h0) / 2; + xe = -b / (2 * a); + ye = (a * xe + b) * xe + h1; + d = b * b - 4 * a * h1; + roots = 0; + + if (d >= 0) { + dx = Math.sqrt(d) / (Math.abs(a) * 2); + x1 = xe - dx; + x2 = xe + dx; + if (Math.abs(x1) <= 1) roots++; + if (Math.abs(x2) <= 1) roots++; + if (x1 < -1) x1 = x2; + } + + if (roots === 1) { + if (h0 < 0) rise = i + x1; + else set = i + x1; + + } else if (roots === 2) { + rise = i + (ye < 0 ? x2 : x1); + set = i + (ye < 0 ? x1 : x2); + } + + if (rise && set) break; + + h0 = h2; + } + + var result = {}; + + if (rise) result.rise = hoursLater(t, rise); + if (set) result.set = hoursLater(t, set); + + if (!rise && !set) result[ye > 0 ? 'alwaysUp' : 'alwaysDown'] = true; + + return result; +}; \ No newline at end of file