-
Notifications
You must be signed in to change notification settings - Fork 115
/
sun-position.js
112 lines (91 loc) · 3.53 KB
/
sun-position.js
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
function sunPosition (year, month, day, hour, minutes, sec, lat, long) {
hour = typeof hour !== 'undefined' ? hour : 12;
minutes = typeof minutes !== 'undefined' ? minutes : 0;
sec = typeof sec !== 'undefined' ? sec : 0;
lat = typeof lat !== 'undefined' ? lat : 46.5;
long = typeof long !== 'undefined' ? long : 6.5;
PI = Math.PI;
pi = Math.PI;
twopi = 2 * pi;
deg2rad = (PI * 2) / 360;
// The input to the Atronomer's almanach is the difference between
// the Julian date and JD 2451545.0 (noon, 1 January 2000)
Date.prototype.getJulian = function() {
return Math.floor((this / 86400000) - (this.getTimezoneOffset()/1440) + 2440587.5);
}
today = new Date(year, month, day, hour, minutes, sec, 0);
var jd = today.getJulian(); //get Julian counterpart
time = jd - 51545;
// console.log(jd);
// Ecliptic coordinates
// Mean longitude
mnlong = 280.460 + 0.9856474 * time;
mnlong = mnlong % 360;
if (mnlong < 0){mnlong = mnlong+360;}
// Mean anomaly
mnanom = 357.528 + 0.9856003 * time;
mnanom = mnanom % 360;
if (mnanom < 0){mnanom = mnanom+360;}
//mnanom[mnanom < 0] = mnanom[mnanom < 0] + 360
mnanom = mnanom * deg2rad
// Ecliptic longitude and obliquity of ecliptic
eclong = mnlong + 1.915 * Math.sin(mnanom) + 0.020 * Math.sin(2 * mnanom);
eclong = eclong % 360;
if (eclong < 0){eclong = eclong+360;}
//eclong[eclong < 0] = eclong[eclong < 0] + 360
oblqec = 23.439 - 0.0000004 * time;
eclong = eclong * deg2rad;
oblqec = oblqec * deg2rad;
// Celestial coordinates
// Right ascension and declination
num = Math.cos(oblqec) * Math.sin(eclong);
den = Math.cos(eclong);
ra = Math.atan(num / den);
if(den<0) {ra = ra+pi;}
//ra[den < 0] = ra[den < 0] + pi
if(den>=0&&num<0) {ra = ra+twopi;}
//ra[den >= 0 & num < 0] = ra[den >= 0 & num < 0] + twopi
dec = Math.asin(Math.sin(oblqec) * Math.sin(eclong));
// Local coordinates
// Greenwich mean sidereal time
gmst = 6.697375 + 0.0657098242 * time + hour;
gmst = gmst % 24;
if (gmst<0) {gmst = gmst +24;}
//gmst[gmst < 0] = gmst[gmst < 0] + 24;
// Local mean sidereal time
lmst = gmst + long / 15;
lmst = lmst % 24;
if (lmst < 0) {lmst = lmst + 24;}
//lmst[lmst < 0] = lmst[lmst < 0] + 24.
lmst = lmst * 15 * deg2rad;
// Hour angle
ha = lmst - ra;
if (ha < -pi) {ha = ha +twopi;}
//ha[ha < -pi] = ha[ha < -pi] + twopi
if (ha > pi) {ha = ha - twopi;}
//ha[ha > pi] = ha[ha > pi] - twopi
// Latitude to radians
lat = lat * deg2rad;
// Azimuth and elevation
el = Math.asin(Math.sin(dec) * Math.sin(lat) + Math.cos(dec) * Math.cos(lat) * Math.cos(ha));
az = Math.asin(-Math.cos(dec) * Math.sin(ha) / Math.cos(el));
// For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353
Math.cosAzPos = (0 <= Math.sin(dec) - Math.sin(el) * Math.sin(lat));
Math.sinAzNeg = (Math.sin(az) < 0);
if (Math.cosAzPos && Math.sinAzNeg) {az = az +twopi;}
//az[Math.cosAzPos & Math.sinAzNeg] = az[Math.cosAzPos & Math.sinAzNeg] + twopi
if (!Math.cosAzPos){az = pi - az;}
//az[!Math.cosAzPos] = pi - az[!Math.cosAzPos]
// if (0 < Math.sin(dec) - Math.sin(el) * Math.sin(lat)) {
// if(Math.sin(az) < 0) {az = az + twopi;}
// } else {
// az = pi - az;
// }
el = el / deg2rad;
az = az / deg2rad;
lat = lat / deg2rad;
//return(list(elevation=el, azimuth=az))
return [el, az];
}
//sunPosition(2012, 06, 23);
// console.log(sunPosition(2012, 06, 23));