計算太陽高度角(javascript實作)

1 /******************************************/
  2 /************** 相關轉換工具 ***************/
  3 /******************************************/
  4 /**
  5  * 
  6  * @param {Number} hour 
  7  * @param {Number} mim 
  8  * @param {Number} sec 
  9  * @param {Boolean} pm 
 10  * @param {Boolean} dst 
 11  * @returns 傳回目前時間過了多少分鐘(0點開始計時)
 12  */
 13 function getLocalTime(hour, mim, sec, pm=false, dst=false) {
 14     if (pm && hour < 12) dochr += 12;
 15     if (dst) dochr -= 1
 16     return hour * 60 + mim + sec / 60.0
 17 }
 18 
 19 // 判斷是否是閏年
 20 function isLeapYear(year) {
 21     return ((year % 4 == 0 && year % 100 != 0) || year % 400 == 0);
 22 }
 23 
 24 // 弧度制轉角度制
 25 function radToDeg(angleRad) {
 26     return 180.0 * angleRad / Math.PI;
 27 }
 28 
 29 // 角度制轉弧度制
 30 function degToRad(angleDeg) {
 31     return Math.PI * angleDeg / 180.0;
 32 }
 33 
 34 // 根據<年,月,日>擷取儒略日
 35 function getJulianDay(year, month, day) {
 36     var A = Math.floor(year / 100)
 37     var B = 2 - A + Math.floor(A / 4)
 38     var JD = Math.floor(365.25 * (year + 4716)) + Math.floor(30.6001 * (month + 1)) + day + B - 1524.5
 39     return JD
 40 }
 41 
 42 /**
 43  * @param {Number} JulianDay => 儒略日
 44  * 儒略日的起始為公元前 4713 年 1 月 1 日中午 12 點
 45  * 每過一“平太陽日”加一 => 即通常意義上的“一天”
 46  * 地球在軌道上做的是不等速運動,一年之内真太陽日的長度不斷改變,于是引進平太陽的概念
 47  * 一年中<真太陽日>的平均即為<平太陽日>
 48  * 2000 年 1 月 1 日的 UT12:00 是儒略日 2451545
 49  * @returns 儒略日至 2000 年始經過的世紀數(一儒略年有 365.25 天)
 50  */
 51 function calcTimeJulianCent(JulianDay) {
 52     return (JulianDay - 2451545.0) / 36525.0;
 53 }
 54 
 55 // 将 d 控制在 0~Range 内
 56 function correctRange(value, Range) {
 57     value = value % Range;
 58     if (value < 0) {
 59         value += Range;
 60     }
 61     return value;
 62 }
 63 
 64 /**
 65  * 将時間(以分鐘為機關)格式化
 66  * @param {Number} minutes 
 67  * @param {Number} flag
 68  * flag = 1 格式化成 HH
 69  * flag = 2 格式化成 HH:MM
 70  * flag = 3 格式化成 HH:MM:SS
 71  */
 72  function timeString(minutes, flag) {
 73     if (minutes >= 0 && minutes < 1440) {
 74         var floatHour = minutes / 60.0;
 75         var hour = Math.floor(floatHour);
 76 
 77         var floatMinute = 60.0 * (floatHour - Math.floor(floatHour));
 78         var minute = Math.floor(floatMinute);
 79 
 80         var floatSec = 60.0 * (floatMinute - Math.floor(floatMinute));
 81         var second = Math.floor(floatSec + 0.5); // 大于 59.5s 就向上取整
 82 
 83         if (second > 59) { // second>59 隻會是 60s
 84             second = 0;
 85             minute += 1;
 86         }
 87         // 需要格式化成 HH:MM 時
 88         if (flag == 2 && second >= 30) minute++; // 當超過 30s 時分鐘數進 1
 89         if (minute > 59) {
 90             minute = 0;
 91             hour += 1;
 92         }
 93         var output = zeroPad(hour, 2) + ':' + zeroPad(minute, 2);
 94         if (flag > 2) output = output + ':' + zeroPad(second, 2);
 95     } else {
 96         var output = 'error'
 97     }
 98     return output;
 99 }
100 
101 // 不足補 0 
102 function zeroPad(n, digits) {
103     n = String(n);
104     while (n.length < digits) {
105         n = '0' + n;
106     }
107     return n;
108 }
109 
110 /******************************************/
111 /********** 太陽坐标計算開始 ***************/
112 /******************************************/
113 
114 /***
115  * @param {Number} t => Julian century(經過的世紀數)
116  * 計算太陽的平黃經(geometric mean longitude [L0])
117  * 黃經：天球黃道坐标系中的經度 => 用來确定天體在天球上的位置
118  */
119 function calcGeomMeanLongSun(t) {
120     let L0 = 280.46646 + 36000.76983 * t + 0.0003032 * Math.pow(t, 2);
121     return correctRange(L0, 360);
122 }
123 
124 /**
125  * 計算太陽的平近點角(mean anomaly [M])
126  * 平近點角：軌道上的物體在輔助圓上相對于中心點的運作角度
127  * 參考：https://en.wikipedia.org/wiki/Mean_anomaly
128  */
129 function calcGeomMeanAnomalySun(t) {
130     let M = 357.52911 + 35999.05029 * t - 0.0001537 * Math.pow(t, 2);
131     return correctRange(M, 360);
132 }
133 
134 /**
135  * 計算地球軌道偏心率(eccentricity [e])
136  * 數學上稱為“離心率”
137  * 對于橢圓 => 兩焦點間的距離(2c)和長軸長度(2a)的比值,即 e=c/a
138  */
139 function calcEccentricityEarthOrbit(t) {
140     let e = 0.016708634 - 0.000042037 * t + 0.0000001267 * Math.pow(t, 2);
141     return e;
142 }
143 
144 /**
145  * 圓心方程(Equation of the Center)
146  * 計算角差(橢圓軌道上實際位置與它在同一周期的圓形軌道上勻速運動時所占位置之間的角差)
147  * 參考 https://en.wikipedia.org/wiki/Equation_of_the_center
148  */
149 function calcSunEqOfCenter(t) {
150     var M = calcGeomMeanAnomalySun(t); // 太陽平近點角
151     var mrad = degToRad(M);
152     var sinm = Math.sin(mrad);
153     var sin2m = Math.sin(2 * mrad);
154     var sin3m = Math.sin(3 * mrad);
155     var C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289;
156     return C; // 機關為度
157 }
158 
159 // 計算太陽真實經度(True Longitude)
160 function calcSunTrueLong(t) {
161     var l0 = calcGeomMeanLongSun(t); // 太陽平黃經
162     var C = calcSunEqOfCenter(t); // 中心方程
163     var Ltrue = correctRange(l0 + C, 360);
164     return Ltrue; // 機關為度
165 }
166 
167 // 計算太陽真實平近點角(True Anomaly)
168 function calcSunTrueAnomaly(t) {
169     var M = calcGeomMeanAnomalySun(t);
170     var C = calcSunEqOfCenter(t);
171     var nu = correctRange(M + C, 360);
172     return nu; // ν(nu)機關為度
173 }
174 
175 /**
176  * 計算半徑矢量 => 太陽中心到地球中心的距離
177  */
178 function calcSunRadVector(t) {
179     var n = calcSunTrueAnomaly(t);
180     var e = calcEccentricityEarthOrbit(t);
181     var R = (1.000001018 * (1 - Math.pow(e, 2))) / (1 + e * Math.cos(degToRad(n)));
182     return R;
183 }
184 
185 /**
186  * 計算太陽表觀經度
187  * 參考 https://en.wikipedia.org/wiki/Apparent_longitude
188  */
189 function calcSunApparentLong(t) {
190     var Ltrue = calcSunTrueLong(t);
191     var omega = 125.04 - 1934.136 * t;
192     var Lapp = Ltrue - 0.00569 - 0.00478 * Math.sin(degToRad(omega));
193     return Lapp; // 機關為度
194 }
195 
196 /**
197  * 計算黃赤交角
198  */
199 function calcMeanObliquityOfEcliptic(t) {
200     var seconds = 21.448 - t * (46.8150 + t * (0.00059 - t * (0.001813)));
201     var e0 = 23.0 + (26.0 + (seconds / 60.0)) / 60.0;
202     return e0; // 機關為度
203 }
204 
205 /**
206  * 計算太陽坐标
207  */
208 function calcObliquityCorrection(t) {
209     var e0 = calcMeanObliquityOfEcliptic(t);
210     var omega = 125.04 - 1934.136 * t;
211     var e = e0 + 0.00256 * Math.cos(degToRad(omega));
212     return e; // 機關為度
213 }
214 
215 // 計算太陽赤經
216 function calcSunRtAscension(t) {
217     var e = calcObliquityCorrection(t);
218     var Lapp = calcSunApparentLong(t);
219     var tananum = (Math.cos(degToRad(e)) * Math.sin(degToRad(Lapp)));
220     var tanadenom = (Math.cos(degToRad(Lapp)));
221     var alpha = radToDeg(Math.atan2(tananum, tanadenom));
222     return alpha.toFixed(2); // 機關為度
223 }
224 
225 // 計算太陽赤緯
226 function calcSunDeclination(t) {
227     var e = calcObliquityCorrection(t);
228     var Lapp = calcSunApparentLong(t);
229     var sint = Math.sin(degToRad(e)) * Math.sin(degToRad(Lapp));
230     var theta = radToDeg(Math.asin(sint));
231     return theta.toFixed(2); // 機關為度
232 }
233 
234 /**
235  * 計算真太陽日與平太陽日之差(即“時差”)
236  * 與目前是幾點無關 => 即計算的是 00:00:00 的時差
237  */
238 function calcEquationOfTime(t) {
239     var epsilon = calcObliquityCorrection(t);
240     var l0 = calcGeomMeanLongSun(t);
241     var e = calcEccentricityEarthOrbit(t);
242     var m = calcGeomMeanAnomalySun(t);
243 
244     var y = Math.tan(degToRad(epsilon) / 2.0);
245     y *= y;
246 
247     var sin2l0 = Math.sin(2.0 * degToRad(l0));
248     var sinm = Math.sin(degToRad(m));
249     var cos2l0 = Math.cos(2.0 * degToRad(l0));
250     var sin4l0 = Math.sin(4.0 * degToRad(l0));
251     var sin2m = Math.sin(2.0 * degToRad(m));
252 
253     var Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m;
254     return Math.floor(radToDeg(Etime) * 4.0 * 100 + 0.5) / 100.0; // 機關為分鐘
255 }
256 
257 /**
258  * 根據目前時間對儒略日進行修正
259  */
260 function CorrectJD(t, hour, mim, sec, tz) {
261     //  擷取分鐘數
262     let mims = getLocalTime(hour, mim, sec);
263     // 擷取時間修正的儒略日(+小時數-時區差)
264     let JD = t + mims / 1440.0 - tz / 24.0;
265     return JD;
266 }
267 
268 // 計算日出角
269 function calcHourAngleSunrise(lat, solarDec) {
270     var latRad = degToRad(lat);
271     var sdRad = degToRad(solarDec);
272     var HAarg = (Math.cos(degToRad(90.833)) / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad) * Math.tan(sdRad));
273     var HA = Math.acos(HAarg);
274     return HA; // 機關弧度
275 }
276 
277 // 計算日落角
278 function calcHourAngleSunset(lat, solarDec) {
279     return -calcHourAngleSunrise(lat, solarDec); // 機關弧度
280 }
281 
282 /******************************************/
283 /********** 太陽坐标計算結束 ***************/
284 /******************************************/
285 
286 /**
287  * 計算地球上給定位置、給定日期“太陽正午”的世界協調時間(UTC)
288  * 太陽高度角最大的時候
289  * @return HH:MM:SS 格式的正午時間
290  */
291 function calcSolNoon(julianday, longitude, timezone, dst) {
292     var tnoon = calcTimeJulianCent(julianday - longitude / 360.0)
293     var eqTime = calcEquationOfTime(tnoon)
294     var solNoonOffset = 720.0 - (longitude * 4) - eqTime // 機關為分鐘
295     var newt = calcTimeJulianCent(jd + solNoonOffset / 1440.0)
296     eqTime = calcEquationOfTime(newt)
297     solNoonLocal = 720 - (longitude * 4) - eqTime + (timezone * 60.0) // 機關為分鐘
298     if (dst) solNoonLocal += 60.0 // 如果采用夏時令則加 60 分鐘(提前 1 小時)
299     // 将時間控制在 24 小時(1440 分鐘) 内
300     solNoonLocal = correctRange(solNoonLocal, 1440.0)
301     return timeString(solNoonLocal, 3);
302 }
303 
304 /**
305  * @param {Number} rise rise = 1 計算日出時間, 0 計算日落時間
306  * @param {Number} JD 儒略日
307  * @param {Number} latitude 緯度
308  * @param {Number} longitude 經度
309  * 傳回日出日落時間的分鐘數
310  */
311 function calcSunriseSetUTC(rise, JD, latitude, longitude) {
312     var t = calcTimeJulianCent(JD);
313     var eqTime = calcEquationOfTime(t);
314     var solarDec = calcSunDeclination(t);
315     var hourAngle = calcHourAngleSunrise(latitude, solarDec);
316     if (!rise) hourAngle = -hourAngle;
317     var delta = longitude + radToDeg(hourAngle);
318     var timeUTC = 720 - (4.0 * delta) - eqTime; // 機關為分鐘
319     return timeUTC;
320 }
321 
322 // 計算日出日落時間
323 function calcSunriseSet(rise, JD, latitude, longitude, timezone, dst) {
324     var timeUTC = calcSunriseSetUTC(rise, JD, latitude, longitude);
325     var newTimeUTC = calcSunriseSetUTC(rise, JD + timeUTC / 1440.0, latitude, longitude);
326     // 時區修正的日出日落分鐘數
327     var timeLocal = newTimeUTC + (timezone * 60.0)
328     // 夏時令修正
329     timeLocal += ((dst) ? 60.0 : 0.0);
330     // 控制 timeLocal 範圍在 0~1440.0 内
331     timeLocal = correctRange(timeLocal, 1440.0);
332     return timeString(timeLocal, 2);   
333 }
334 
335 /**
336  * @param {Number} localtime => 通過 getLocalTime 獲得的時間
337  * 計算太陽高度角
338  */
339 function calcelevation(t, localtime, latitude, longitude, zone) {
340     var eqTime = calcEquationOfTime(t); // 時差
341     var theta = calcSunDeclination(t); // 赤緯
342     var solarTimeFix = eqTime + 4.0 * longitude - 60.0 * zone;
343 
344     var trueSolarTime = localtime + solarTimeFix;
345     
346     trueSolarTime = correctRange(trueSolarTime, 1440);
347 
348     var hourAngle = trueSolarTime / 4.0 - 180.0;
349     if (hourAngle < -180) hourAngle += 360.0;
350 
351     var haRad = degToRad(hourAngle);
352     var csz = Math.sin(degToRad(latitude)) * Math.sin(degToRad(theta)) + Math.cos(degToRad(latitude)) * Math.cos(degToRad(theta)) * Math.cos(haRad)
353     if (csz > 1.0) {
354         csz = 1.0
355     } else if (csz < -1.0) {
356         csz = -1.0
357     }
358     var zenith = radToDeg(Math.acos(csz))
359 
360     var exoatmElevation = 90.0 - zenith
361 
362     if (exoatmElevation > 85.0) {
363         var refractionCorrection = 0.0;
364     } else {
365         var te = Math.tan(degToRad(exoatmElevation));
366         if (exoatmElevation > 5.0) {
367             var refractionCorrection = 58.1 / te - 0.07 / (te * te * te) + 0.000086 / (te * te * te * te * te);
368         } else if (exoatmElevation > -0.575) {
369             var refractionCorrection = 1735.0 + exoatmElevation * (-518.2 + exoatmElevation * (103.4 + exoatmElevation * (-12.79 + exoatmElevation * 0.711)));
370         } else {
371             var refractionCorrection = -20.774 / te;
372         }
373         refractionCorrection = refractionCorrection / 3600.0;
374     }
375 
376     var solarZen = zenith - refractionCorrection;
377 
378     if (solarZen > 108.0) {
379         document.getElementById("azbox").value = "dark";
380         document.getElementById("elbox").value = "dark";
381     } else {
382         // 太陽高度角
383         return Math.floor((90.0 - solarZen) * 100 + 0.5) / 100.0
384     }
385 }
386 
387 /**
388  * @param {Number} ZeroAzimuth 零方位角
389  * 零方位角 = 北, ZeroAzimuth=0
390  * 零方位角 = 南, ZeroAzimuth=180
391  */
392 function calcazimuth(t, localtime, latitude, longitude, zone, ZeroAzimuth) {
393     var eqTime = calcEquationOfTime(t); // 時差
394     var theta = calcSunDeclination(t); // 赤緯
395     var solarTimeFix = eqTime + 4.0 * longitude - 60.0 * zone;
396 
397     var trueSolarTime = localtime + solarTimeFix;
398     
399     trueSolarTime = correctRange(trueSolarTime, 1440);
400 
401     var hourAngle = trueSolarTime / 4.0 - 180.0;
402     if (hourAngle < -180) hourAngle += 360.0;
403 
404     var haRad = degToRad(hourAngle);
405     var csz = Math.sin(degToRad(latitude)) * Math.sin(degToRad(theta)) + Math.cos(degToRad(latitude)) * Math.cos(degToRad(theta)) * Math.cos(haRad)
406     if (csz > 1.0) {
407         csz = 1.0
408     } else if (csz < -1.0) {
409         csz = -1.0
410     }
411     var zenith = radToDeg(Math.acos(csz))
412     var azDenom = (Math.cos(degToRad(latitude)) * Math.sin(degToRad(zenith)))
413     
414     if (Math.abs(azDenom) > 0.001) {
415         azRad = ((Math.sin(degToRad(latitude)) * Math.cos(degToRad(zenith))) - Math.sin(degToRad(theta))) / azDenom
416         if (Math.abs(azRad) > 1.0) {
417             if (azRad < 0) {
418                 azRad = -1.0
419             } else {
420                 azRad = 1.0
421             }
422         }
423         var azimuth = 180.0 - radToDeg(Math.acos(azRad))
424         if (hourAngle > 0.0) {
425             azimuth = -azimuth
426         }
427     } else {
428         if (latitude > 0.0) {
429             azimuth = 180.0
430         } else {
431             azimuth = 0.0
432         }
433     }
434     if (azimuth < 0.0) {
435         azimuth += 360.0
436     }
437     return (Math.floor(azimuth * 100 + 0.5) - ZeroAzimuth * 100) / 100.0;
438 }      

相關方法      

• getLocalTime　　計算今天過了多少分鐘
• calcSunriseSet 　　計算日出日落時間
• calcSolNoon 　　計算正午時間
• calcSunRtAscension　　計算太陽赤經(未進行時間修正)
• calcSunDeclination　　計算太陽赤緯(未進行時間修正)
• calcSunRadVector　　計算矢量半徑(未進行時間修正)
• calcMeanObliquityOfEcliptic　　計算黃赤交角(未進行時間修正)
• calcEquationOfTime　　計算時差(未進行時間修正)
• getJulianDay　　擷取儒略日
• CorrectJD　　時間修正的儒略日
• calcTimeJulianCent　　根據儒略日擷取世紀數
• calcelevation　　計算太陽高度角
• calcazimuth　　計算方位角

使用

　

// 擷取目前時間
let date = new Date()
let year = date.getFullYear()
let month = date.getMonth()
let day = date.getDate()

let hour = date.getHours()
let mim = date.getMinutes()
let second = date.getSeconds()

// 擷取儒略日
let jd = getJulianDay(year, month+1, day);

// 經緯度
let latitude = 30.214512
let longitude = 120.225365

// 時區
let tz = 8

// 今天已過了多少分鐘
let tl = getLocalTime(hour,mim,second)

// 使用時間、時區修正的分鐘數
let t = jd + tl / 1440.0 - tz / 24.0

// 世紀數
t =  calcTimeJulianCent(t)

console.log(calcazimuth(t,tl, latitude, longitude,tz, 0)) // 方位角
console.log(calcelevation(t,tl, latitude, longitude,tz)) // 高度角
console.log(calcEquationOfTime(t)) // 時差
console.log(calcSunRtAscension(t)) // 赤經
console.log(calcSunDeclination(t)) //赤緯      

驗證位址  https://gml.noaa.gov/grad/solcalc/

​​ESRL Global Monitoring Laboratory - Global Radiation and Aerosols (noaa.gov)​​