計算太陽高度角(javascript實現)
阿新 • • 發佈：2021-08-30
程式碼來源 https://gitee.com/zhonghongqiang/calculateSunAzEl.git
程式碼參考 https://www.jianshu.com/p/ecdccbfe0a88
太陽座標計算參考資料 https://squarewidget.com/solar-coordinates/
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; 15if (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 } 2829 // 角度制轉弧度制 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)