Browse Source

add more accurate sunset/rise times

Kyle Perik 1 year ago
parent
commit
eeff91f55b
3 changed files with 342 additions and 7 deletions
  1. 1 0
      index.html
  2. 31 7
      index.js
  3. 310 0
      suncalc.js

+ 1 - 0
index.html

@@ -4,6 +4,7 @@
4 4
     </head>
5 5
     <body style="margin: 0; overflow: hidden;">
6 6
         <canvas id="pixelCanvas" width="800" height="800"></canvas>
7
+        <script src="suncalc.js"></script>
7 8
         <script src="index.js"></script>
8 9
     </body>
9 10
 </html>

+ 31 - 7
index.js

@@ -9,6 +9,27 @@ function getPixelCanvas(ctx, size) {
9 9
     };
10 10
 }
11 11
 
12
+function getSunPosition(coords) {
13
+    return SunCalc.getPosition(new Date(), coords.lat, coords.long).altitude;
14
+}
15
+
16
+function getLatLong() {
17
+    return new Promise((resolve, reject) => {
18
+        //navigator.geolocation.getCurrentPosition(resolve, reject);
19
+        resolve({ coords: { latitude: 43, longitude: - 88 } });
20
+    }).then(
21
+        r => ({ lat: r.coords.latitude, long: r.coords.longitude })
22
+    ).catch(
23
+        r => ({ lat: 43, long: -88 })
24
+    );
25
+}
26
+
27
+function getCurrentSunPos() {
28
+    return getLatLong().then(
29
+        r => getSunPosition(r)
30
+    );
31
+}
32
+
12 33
 function drawRect(x, y, width, height, color, pc) {
13 34
     pc.ctx.fillStyle = color;
14 35
     pc.ctx.fillRect(x * pc.size, y * pc.size, width * pc.size, height * pc.size);
@@ -43,11 +64,11 @@ var starcolor = [255, 255, 255];
43 64
 var sunsetColor = [255, 100, 50];
44 65
 var sunset2Color = [255, 50, 50];
45 66
 
46
-function getcolorat(x, y, i, width, height, isstar) {
67
+function getcolorat(x, y, percentday, width, height, isstar) {
47 68
     var star = isstar ? starcolor : null;
48 69
     var today = new Date();
49 70
     var daypercentcomplete = (today.getHours() + today.getMinutes() / 60 + today.getSeconds() / 3600) / 24;
50
-    var percentday = Math.sin(daypercentcomplete * Math.PI * 2 - Math.PI / 2);
71
+    // var percentday = Math.sin(daypercentcomplete * Math.PI * 2 - Math.PI / 2);
51 72
     var sunheight = percentday * height;
52 73
     var sunpos = [width / 2, sunheight - height];
53 74
     var sunshineColor = colorMix(
@@ -92,7 +113,7 @@ function getSize(pc) {
92 113
 }
93 114
 
94 115
 function draw(s) {
95
-    var i = s.i;
116
+    var percentday = getSunPosition(s.coords);
96 117
     var canvasSize = getSize(pixelCanvas);
97 118
     drawRect(0, 0, canvasSize.x, canvasSize.y, 'white', pixelCanvas);
98 119
     [...Array(canvasSize.x).keys()].map(
@@ -100,13 +121,13 @@ function draw(s) {
100 121
             x, y,
101 122
             1, 1,
102 123
             rgbtohex(getcolorat(
103
-            	x, y, i, canvasSize.x, canvasSize.y,
124
+            	x, y, percentday, canvasSize.x, canvasSize.y,
104 125
             	s.stars.includes(x + y * canvasSize.x)
105 126
             )),
106 127
             pixelCanvas
107 128
         ))
108 129
     );
109
-    return {i: i + 1, stars: s.stars};
130
+    return s;
110 131
 }
111 132
 
112 133
 
@@ -123,7 +144,7 @@ function blowup() {
123 144
 
124 145
 var size = getSize(pixelCanvas);
125 146
 var state = {
126
-    i: 0,
147
+    coords: {},
127 148
     stars: [...Array(100).keys()].map(
128 149
 	n =>  Math.round(
129 150
 	    Math.random() * size.x * size.y
@@ -131,6 +152,8 @@ var state = {
131 152
     )
132 153
 };
133 154
 
155
+getLatLong().then(r => state.coords = r);
156
+
134 157
 setInterval(() => state = draw(state), 1000);
135 158
 
136 159
 blowup();
@@ -138,10 +161,11 @@ blowup();
138 161
 var resizeTimeout;
139 162
 window.addEventListener('resize', () => {
140 163
     clearTimeout(resizeTimeout);
141
-    resizeTimeout = setTimeout(blowup, 100);
164
+    resizeTimeout = setTimeout(blowup, 10000);
142 165
     state.stars = [...Array(100).keys()].map(
143 166
         n => Math.round(
144 167
             Math.random() * canvas.width * canvas.height
145 168
         )
146 169
     );
170
+    draw(state);
147 171
 });

+ 310 - 0
suncalc.js

@@ -0,0 +1,310 @@
1
+/*
2
+ (c) 2011-2015, Vladimir Agafonkin
3
+ SunCalc is a JavaScript library for calculating sun/moon position and light phases.
4
+ https://github.com/mourner/suncalc
5
+*/
6
+
7
+(function () { 'use strict';
8
+
9
+// shortcuts for easier to read formulas
10
+
11
+var PI   = Math.PI,
12
+    sin  = Math.sin,
13
+    cos  = Math.cos,
14
+    tan  = Math.tan,
15
+    asin = Math.asin,
16
+    atan = Math.atan2,
17
+    acos = Math.acos,
18
+    rad  = PI / 180;
19
+
20
+// sun calculations are based on http://aa.quae.nl/en/reken/zonpositie.html formulas
21
+
22
+
23
+// date/time constants and conversions
24
+
25
+var dayMs = 1000 * 60 * 60 * 24,
26
+    J1970 = 2440588,
27
+    J2000 = 2451545;
28
+
29
+function toJulian(date) { return date.valueOf() / dayMs - 0.5 + J1970; }
30
+function fromJulian(j)  { return new Date((j + 0.5 - J1970) * dayMs); }
31
+function toDays(date)   { return toJulian(date) - J2000; }
32
+
33
+
34
+// general calculations for position
35
+
36
+var e = rad * 23.4397; // obliquity of the Earth
37
+
38
+function rightAscension(l, b) { return atan(sin(l) * cos(e) - tan(b) * sin(e), cos(l)); }
39
+function declination(l, b)    { return asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l)); }
40
+
41
+function azimuth(H, phi, dec)  { return atan(sin(H), cos(H) * sin(phi) - tan(dec) * cos(phi)); }
42
+function altitude(H, phi, dec) { return asin(sin(phi) * sin(dec) + cos(phi) * cos(dec) * cos(H)); }
43
+
44
+function siderealTime(d, lw) { return rad * (280.16 + 360.9856235 * d) - lw; }
45
+
46
+function astroRefraction(h) {
47
+    if (h < 0) // the following formula works for positive altitudes only.
48
+        h = 0; // if h = -0.08901179 a div/0 would occur.
49
+
50
+    // formula 16.4 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
51
+    // 1.02 / tan(h + 10.26 / (h + 5.10)) h in degrees, result in arc minutes -> converted to rad:
52
+    return 0.0002967 / Math.tan(h + 0.00312536 / (h + 0.08901179));
53
+}
54
+
55
+// general sun calculations
56
+
57
+function solarMeanAnomaly(d) { return rad * (357.5291 + 0.98560028 * d); }
58
+
59
+function eclipticLongitude(M) {
60
+
61
+    var C = rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M)), // equation of center
62
+        P = rad * 102.9372; // perihelion of the Earth
63
+
64
+    return M + C + P + PI;
65
+}
66
+
67
+function sunCoords(d) {
68
+
69
+    var M = solarMeanAnomaly(d),
70
+        L = eclipticLongitude(M);
71
+
72
+    return {
73
+        dec: declination(L, 0),
74
+        ra: rightAscension(L, 0)
75
+    };
76
+}
77
+
78
+
79
+var SunCalc = {};
80
+
81
+
82
+// calculates sun position for a given date and latitude/longitude
83
+
84
+SunCalc.getPosition = function (date, lat, lng) {
85
+
86
+    var lw  = rad * -lng,
87
+        phi = rad * lat,
88
+        d   = toDays(date),
89
+
90
+        c  = sunCoords(d),
91
+        H  = siderealTime(d, lw) - c.ra;
92
+
93
+    return {
94
+        azimuth: azimuth(H, phi, c.dec),
95
+        altitude: altitude(H, phi, c.dec)
96
+    };
97
+};
98
+
99
+
100
+// sun times configuration (angle, morning name, evening name)
101
+
102
+var times = SunCalc.times = [
103
+    [-0.833, 'sunrise',       'sunset'      ],
104
+    [  -0.3, 'sunriseEnd',    'sunsetStart' ],
105
+    [    -6, 'dawn',          'dusk'        ],
106
+    [   -12, 'nauticalDawn',  'nauticalDusk'],
107
+    [   -18, 'nightEnd',      'night'       ],
108
+    [     6, 'goldenHourEnd', 'goldenHour'  ]
109
+];
110
+
111
+// adds a custom time to the times config
112
+
113
+SunCalc.addTime = function (angle, riseName, setName) {
114
+    times.push([angle, riseName, setName]);
115
+};
116
+
117
+
118
+// calculations for sun times
119
+
120
+var J0 = 0.0009;
121
+
122
+function julianCycle(d, lw) { return Math.round(d - J0 - lw / (2 * PI)); }
123
+
124
+function approxTransit(Ht, lw, n) { return J0 + (Ht + lw) / (2 * PI) + n; }
125
+function solarTransitJ(ds, M, L)  { return J2000 + ds + 0.0053 * sin(M) - 0.0069 * sin(2 * L); }
126
+
127
+function hourAngle(h, phi, d) { return acos((sin(h) - sin(phi) * sin(d)) / (cos(phi) * cos(d))); }
128
+
129
+// returns set time for the given sun altitude
130
+function getSetJ(h, lw, phi, dec, n, M, L) {
131
+
132
+    var w = hourAngle(h, phi, dec),
133
+        a = approxTransit(w, lw, n);
134
+    return solarTransitJ(a, M, L);
135
+}
136
+
137
+
138
+// calculates sun times for a given date and latitude/longitude
139
+
140
+SunCalc.getTimes = function (date, lat, lng) {
141
+
142
+    var lw = rad * -lng,
143
+        phi = rad * lat,
144
+
145
+        d = toDays(date),
146
+        n = julianCycle(d, lw),
147
+        ds = approxTransit(0, lw, n),
148
+
149
+        M = solarMeanAnomaly(ds),
150
+        L = eclipticLongitude(M),
151
+        dec = declination(L, 0),
152
+
153
+        Jnoon = solarTransitJ(ds, M, L),
154
+
155
+        i, len, time, Jset, Jrise;
156
+
157
+
158
+    var result = {
159
+        solarNoon: fromJulian(Jnoon),
160
+        nadir: fromJulian(Jnoon - 0.5)
161
+    };
162
+
163
+    for (i = 0, len = times.length; i < len; i += 1) {
164
+        time = times[i];
165
+
166
+        Jset = getSetJ(time[0] * rad, lw, phi, dec, n, M, L);
167
+        Jrise = Jnoon - (Jset - Jnoon);
168
+
169
+        result[time[1]] = fromJulian(Jrise);
170
+        result[time[2]] = fromJulian(Jset);
171
+    }
172
+
173
+    return result;
174
+};
175
+
176
+
177
+// moon calculations, based on http://aa.quae.nl/en/reken/hemelpositie.html formulas
178
+
179
+function moonCoords(d) { // geocentric ecliptic coordinates of the moon
180
+
181
+    var L = rad * (218.316 + 13.176396 * d), // ecliptic longitude
182
+        M = rad * (134.963 + 13.064993 * d), // mean anomaly
183
+        F = rad * (93.272 + 13.229350 * d),  // mean distance
184
+
185
+        l  = L + rad * 6.289 * sin(M), // longitude
186
+        b  = rad * 5.128 * sin(F),     // latitude
187
+        dt = 385001 - 20905 * cos(M);  // distance to the moon in km
188
+
189
+    return {
190
+        ra: rightAscension(l, b),
191
+        dec: declination(l, b),
192
+        dist: dt
193
+    };
194
+}
195
+
196
+SunCalc.getMoonPosition = function (date, lat, lng) {
197
+
198
+    var lw  = rad * -lng,
199
+        phi = rad * lat,
200
+        d   = toDays(date),
201
+
202
+        c = moonCoords(d),
203
+        H = siderealTime(d, lw) - c.ra,
204
+        h = altitude(H, phi, c.dec),
205
+        // formula 14.1 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
206
+        pa = atan(sin(H), tan(phi) * cos(c.dec) - sin(c.dec) * cos(H));
207
+
208
+    h = h + astroRefraction(h); // altitude correction for refraction
209
+
210
+    return {
211
+        azimuth: azimuth(H, phi, c.dec),
212
+        altitude: h,
213
+        distance: c.dist,
214
+        parallacticAngle: pa
215
+    };
216
+};
217
+
218
+
219
+// calculations for illumination parameters of the moon,
220
+// based on http://idlastro.gsfc.nasa.gov/ftp/pro/astro/mphase.pro formulas and
221
+// Chapter 48 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
222
+
223
+SunCalc.getMoonIllumination = function (date) {
224
+
225
+    var d = toDays(date || new Date()),
226
+        s = sunCoords(d),
227
+        m = moonCoords(d),
228
+
229
+        sdist = 149598000, // distance from Earth to Sun in km
230
+
231
+        phi = acos(sin(s.dec) * sin(m.dec) + cos(s.dec) * cos(m.dec) * cos(s.ra - m.ra)),
232
+        inc = atan(sdist * sin(phi), m.dist - sdist * cos(phi)),
233
+        angle = atan(cos(s.dec) * sin(s.ra - m.ra), sin(s.dec) * cos(m.dec) -
234
+                cos(s.dec) * sin(m.dec) * cos(s.ra - m.ra));
235
+
236
+    return {
237
+        fraction: (1 + cos(inc)) / 2,
238
+        phase: 0.5 + 0.5 * inc * (angle < 0 ? -1 : 1) / Math.PI,
239
+        angle: angle
240
+    };
241
+};
242
+
243
+
244
+function hoursLater(date, h) {
245
+    return new Date(date.valueOf() + h * dayMs / 24);
246
+}
247
+
248
+// calculations for moon rise/set times are based on http://www.stargazing.net/kepler/moonrise.html article
249
+
250
+SunCalc.getMoonTimes = function (date, lat, lng, inUTC) {
251
+    var t = new Date(date);
252
+    if (inUTC) t.setUTCHours(0, 0, 0, 0);
253
+    else t.setHours(0, 0, 0, 0);
254
+
255
+    var hc = 0.133 * rad,
256
+        h0 = SunCalc.getMoonPosition(t, lat, lng).altitude - hc,
257
+        h1, h2, rise, set, a, b, xe, ye, d, roots, x1, x2, dx;
258
+
259
+    // go in 2-hour chunks, each time seeing if a 3-point quadratic curve crosses zero (which means rise or set)
260
+    for (var i = 1; i <= 24; i += 2) {
261
+        h1 = SunCalc.getMoonPosition(hoursLater(t, i), lat, lng).altitude - hc;
262
+        h2 = SunCalc.getMoonPosition(hoursLater(t, i + 1), lat, lng).altitude - hc;
263
+
264
+        a = (h0 + h2) / 2 - h1;
265
+        b = (h2 - h0) / 2;
266
+        xe = -b / (2 * a);
267
+        ye = (a * xe + b) * xe + h1;
268
+        d = b * b - 4 * a * h1;
269
+        roots = 0;
270
+
271
+        if (d >= 0) {
272
+            dx = Math.sqrt(d) / (Math.abs(a) * 2);
273
+            x1 = xe - dx;
274
+            x2 = xe + dx;
275
+            if (Math.abs(x1) <= 1) roots++;
276
+            if (Math.abs(x2) <= 1) roots++;
277
+            if (x1 < -1) x1 = x2;
278
+        }
279
+
280
+        if (roots === 1) {
281
+            if (h0 < 0) rise = i + x1;
282
+            else set = i + x1;
283
+
284
+        } else if (roots === 2) {
285
+            rise = i + (ye < 0 ? x2 : x1);
286
+            set = i + (ye < 0 ? x1 : x2);
287
+        }
288
+
289
+        if (rise && set) break;
290
+
291
+        h0 = h2;
292
+    }
293
+
294
+    var result = {};
295
+
296
+    if (rise) result.rise = hoursLater(t, rise);
297
+    if (set) result.set = hoursLater(t, set);
298
+
299
+    if (!rise && !set) result[ye > 0 ? 'alwaysUp' : 'alwaysDown'] = true;
300
+
301
+    return result;
302
+};
303
+
304
+
305
+// export as Node module / AMD module / browser variable
306
+if (typeof exports === 'object' && typeof module !== 'undefined') module.exports = SunCalc;
307
+else if (typeof define === 'function' && define.amd) define(SunCalc);
308
+else window.SunCalc = SunCalc;
309
+
310
+}());