计算太阳高度角(javascript实现)

 

代码来源 https://gitee.com/zhonghongqiang/calculateSunAzEl.git
代码参考 https://www.jianshu.com/p/ecdccbfe0a88
太阳坐标计算参考资料 https://squarewidget.com/solar-coordinates/

/******************************************/
/************** 相关转换工具 ***************/
/******************************************/
/**
 * 
 * @param {Number} hour 
 * @param {Number} mim 
 * @param {Number} sec 
 * @param {Boolean} pm 
 * @param {Boolean} dst 
 * @returns 返回当前时间过了多少分钟(0点开始计时)
 */
function getLocalTime(hour, mim, sec, pm=false, dst=false) {
    if (pm && hour < 12) dochr += 12;
    if (dst) dochr -= 1
    return hour * 60 + mim + sec / 60.0
}

// 判断是否是闰年
function isLeapYear(year) {
    return ((year % 4 == 0 && year % 100 != 0) || year % 400 == 0);
}

// 弧度制转角度制
function radToDeg(angleRad) {
    return 180.0 * angleRad / Math.PI;
}

// 角度制转弧度制
function degToRad(angleDeg) {
    return Math.PI * angleDeg / 180.0;
}

// 根据<年,月,日>获取儒略日
function getJulianDay(year, month, day) {
    var A = Math.floor(year / 100)
    var B = 2 - A + Math.floor(A / 4)
    var JD = Math.floor(365.25 * (year + 4716)) + Math.floor(30.6001 * (month + 1)) + day + B - 1524.5
    return JD
}

/**
 * @param {Number} JulianDay => 儒略日
 * 儒略日的起始为公元前 4713 年 1 月 1 日中午 12 点
 * 每过一“平太阳日”加一 => 即通常意义上的“一天”
 * 地球在轨道上做的是不等速运动,一年之内真太阳日的长度不断改变,于是引进平太阳的概念
 * 一年中<真太阳日>的平均即为<平太阳日>
 * 2000 年 1 月 1 日的 UT12:00 是儒略日 2451545
 * @returns 儒略日至 2000 年始经过的世纪数(一儒略年有 365.25 天)
 */
function calcTimeJulianCent(JulianDay) {
    return (JulianDay - 2451545.0) / 36525.0;
}

// 将 d 控制在 0~Range 内
function correctRange(value, Range) {
    value = value % Range;
    if (value < 0) {
        value += Range;
    }
    return value;
}

/**
 * 将时间(以分钟为单位)格式化
 * @param {Number} minutes 
 * @param {Number} flag
 * flag = 1 格式化成 HH
 * flag = 2 格式化成 HH:MM
 * flag = 3 格式化成 HH:MM:SS
 */
 function timeString(minutes, flag) {
    if (minutes >= 0 && minutes < 1440) {
        var floatHour = minutes / 60.0;
        var hour = Math.floor(floatHour);

        var floatMinute = 60.0 * (floatHour - Math.floor(floatHour));
        var minute = Math.floor(floatMinute);

        var floatSec = 60.0 * (floatMinute - Math.floor(floatMinute));
        var second = Math.floor(floatSec + 0.5); // 大于 59.5s 就向上取整

        if (second > 59) { // second>59 只会是 60s
            second = 0;
            minute += 1;
        }
        // 需要格式化成 HH:MM 时
        if (flag == 2 && second >= 30) minute++; // 当超过 30s 时分钟数进 1
        if (minute > 59) {
            minute = 0;
            hour += 1;
        }
        var output = zeroPad(hour, 2) + ':' + zeroPad(minute, 2);
        if (flag > 2) output = output + ':' + zeroPad(second, 2);
    } else {
        var output = 'error'
    }
    return output;
}

// 不足补 0 
function zeroPad(n, digits) {
    n = String(n);
    while (n.length < digits) {
        n = '0' + n;
    }
    return n;
}

/******************************************/
/********** 太阳坐标计算开始 ***************/
/******************************************/

/***
 * @param {Number} t => Julian century(经过的世纪数)
 * 计算太阳的平黄经(geometric mean longitude [L0])
 * 黄经：天球黄道坐标系中的经度 => 用来确定天体在天球上的位置
 */
function calcGeomMeanLongSun(t) {
    let L0 = 280.46646 + 36000.76983 * t + 0.0003032 * Math.pow(t, 2);
    return correctRange(L0, 360);
}

/**
 * 计算太阳的平近点角(mean anomaly [M])
 * 平近点角：轨道上的物体在辅助圆上相对于中心点的运行角度
 * 参考：https://en.wikipedia.org/wiki/Mean_anomaly
 */
function calcGeomMeanAnomalySun(t) {
    let M = 357.52911 + 35999.05029 * t - 0.0001537 * Math.pow(t, 2);
    return correctRange(M, 360);
}

/**
 * 计算地球轨道偏心率(eccentricity [e])
 * 数学上称为“离心率”
 * 对于椭圆 => 两焦点间的距离(2c)和长轴长度(2a)的比值,即 e=c/a
 */
function calcEccentricityEarthOrbit(t) {
    let e = 0.016708634 - 0.000042037 * t + 0.0000001267 * Math.pow(t, 2);
    return e;
}

/**
 * 圆心方程(Equation of the Center)
 * 计算角差(椭圆轨道上实际位置与它在同一周期的圆形轨道上匀速运动时所占位置之间的角差)
 * 参考 https://en.wikipedia.org/wiki/Equation_of_the_center
 */
function calcSunEqOfCenter(t) {
    var M = calcGeomMeanAnomalySun(t); // 太阳平近点角
    var mrad = degToRad(M);
    var sinm = Math.sin(mrad);
    var sin2m = Math.sin(2 * mrad);
    var sin3m = Math.sin(3 * mrad);
    var C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289;
    return C; // 单位为度
}

// 计算太阳真实经度(True Longitude)
function calcSunTrueLong(t) {
    var l0 = calcGeomMeanLongSun(t); // 太阳平黄经
    var C = calcSunEqOfCenter(t); // 中心方程
    var Ltrue = correctRange(l0 + C, 360);
    return Ltrue; // 单位为度
}

// 计算太阳真实平近点角(True Anomaly)
function calcSunTrueAnomaly(t) {
    var M = calcGeomMeanAnomalySun(t);
    var C = calcSunEqOfCenter(t);
    var nu = correctRange(M + C, 360);
    return nu; // ν(nu)单位为度
}

/**
 * 计算半径矢量 => 太阳中心到地球中心的距离
 */
function calcSunRadVector(t) {
    var n = calcSunTrueAnomaly(t);
    var e = calcEccentricityEarthOrbit(t);
    var R = (1.000001018 * (1 - Math.pow(e, 2))) / (1 + e * Math.cos(degToRad(n)));
    return R;
}

/**
 * 计算太阳表观经度
 * 参考 https://en.wikipedia.org/wiki/Apparent_longitude
 */
function calcSunApparentLong(t) {
    var Ltrue = calcSunTrueLong(t);
    var omega = 125.04 - 1934.136 * t;
    var Lapp = Ltrue - 0.00569 - 0.00478 * Math.sin(degToRad(omega));
    return Lapp; // 单位为度
}

/**
 * 计算黄赤交角
 */
function calcMeanObliquityOfEcliptic(t) {
    var seconds = 21.448 - t * (46.8150 + t * (0.00059 - t * (0.001813)));
    var e0 = 23.0 + (26.0 + (seconds / 60.0)) / 60.0;
    return e0; // 单位为度
}

/**
 * 计算太阳坐标
 */
function calcObliquityCorrection(t) {
    var e0 = calcMeanObliquityOfEcliptic(t);
    var omega = 125.04 - 1934.136 * t;
    var e = e0 + 0.00256 * Math.cos(degToRad(omega));
    return e; // 单位为度
}

// 计算太阳赤经
function calcSunRtAscension(t) {
    var e = calcObliquityCorrection(t);
    var Lapp = calcSunApparentLong(t);
    var tananum = (Math.cos(degToRad(e)) * Math.sin(degToRad(Lapp)));
    var tanadenom = (Math.cos(degToRad(Lapp)));
    var alpha = radToDeg(Math.atan2(tananum, tanadenom));
    return alpha.toFixed(2); // 单位为度
}

// 计算太阳赤纬
function calcSunDeclination(t) {
    var e = calcObliquityCorrection(t);
    var Lapp = calcSunApparentLong(t);
    var sint = Math.sin(degToRad(e)) * Math.sin(degToRad(Lapp));
    var theta = radToDeg(Math.asin(sint));
    return theta.toFixed(2); // 单位为度
}

/**
 * 计算真太阳日与平太阳日之差(即“时差”)
 * 与当前是几点无关 => 即计算的是 00:00:00 的时差
 */
function calcEquationOfTime(t) {
    var epsilon = calcObliquityCorrection(t);
    var l0 = calcGeomMeanLongSun(t);
    var e = calcEccentricityEarthOrbit(t);
    var m = calcGeomMeanAnomalySun(t);

    var y = Math.tan(degToRad(epsilon) / 2.0);
    y *= y;

    var sin2l0 = Math.sin(2.0 * degToRad(l0));
    var sinm = Math.sin(degToRad(m));
    var cos2l0 = Math.cos(2.0 * degToRad(l0));
    var sin4l0 = Math.sin(4.0 * degToRad(l0));
    var sin2m = Math.sin(2.0 * degToRad(m));

    var Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m;
    return Math.floor(radToDeg(Etime) * 4.0 * 100 + 0.5) / 100.0; // 单位为分钟
}

/**
 * 根据当前时间对儒略日进行修正
 */
function CorrectJD(t, hour, mim, sec, tz) {
    //  获取分钟数
    let mims = getLocalTime(hour, mim, sec);
    // 获取时间修正的儒略日(+小时数-时区差)
    let JD = t + mims / 1440.0 - tz / 24.0;
    return JD;
}

// 计算日出角
function calcHourAngleSunrise(lat, solarDec) {
    var latRad = degToRad(lat);
    var sdRad = degToRad(solarDec);
    var HAarg = (Math.cos(degToRad(90.833)) / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad) * Math.tan(sdRad));
    var HA = Math.acos(HAarg);
    return HA; // 单位弧度
}

// 计算日落角
function calcHourAngleSunset(lat, solarDec) {
    return -calcHourAngleSunrise(lat, solarDec); // 单位弧度
}

/******************************************/
/********** 太阳坐标计算结束 ***************/
/******************************************/

/**
 * 计算地球上给定位置、给定日期“太阳正午”的世界协调时间(UTC)
 * 太阳高度角最大的时候
 * @return HH:MM:SS 格式的正午时间
 */
function calcSolNoon(julianday, longitude, timezone, dst) {
    var tnoon = calcTimeJulianCent(julianday - longitude / 360.0)
    var eqTime = calcEquationOfTime(tnoon)
    var solNoonOffset = 720.0 - (longitude * 4) - eqTime // 单位为分钟
    var newt = calcTimeJulianCent(jd + solNoonOffset / 1440.0)
    eqTime = calcEquationOfTime(newt)
    solNoonLocal = 720 - (longitude * 4) - eqTime + (timezone * 60.0) // 单位为分钟
    if (dst) solNoonLocal += 60.0 // 如果采用夏时令则加 60 分钟(提前 1 小时)
    // 将时间控制在 24 小时(1440 分钟) 内
    solNoonLocal = correctRange(solNoonLocal, 1440.0)
    return timeString(solNoonLocal, 3);
}

/**
 * @param {Number} rise rise = 1 计算日出时间, 0 计算日落时间
 * @param {Number} JD 儒略日
 * @param {Number} latitude 纬度
 * @param {Number} longitude 经度
 * 返回日出日落时间的分钟数
 */
function calcSunriseSetUTC(rise, JD, latitude, longitude) {
    var t = calcTimeJulianCent(JD);
    var eqTime = calcEquationOfTime(t);
    var solarDec = calcSunDeclination(t);
    var hourAngle = calcHourAngleSunrise(latitude, solarDec);
    if (!rise) hourAngle = -hourAngle;
    var delta = longitude + radToDeg(hourAngle);
    var timeUTC = 720 - (4.0 * delta) - eqTime; // 单位为分钟
    return timeUTC;
}

// 计算日出日落时间
function calcSunriseSet(rise, JD, latitude, longitude, timezone, dst) {
    var timeUTC = calcSunriseSetUTC(rise, JD, latitude, longitude);
    var newTimeUTC = calcSunriseSetUTC(rise, JD + timeUTC / 1440.0, latitude, longitude);
    // 时区修正的日出日落分钟数
    var timeLocal = newTimeUTC + (timezone * 60.0)
    // 夏时令修正
    timeLocal += ((dst) ? 60.0 : 0.0);
    // 控制 timeLocal 范围在 0~1440.0 内
    timeLocal = correctRange(timeLocal, 1440.0);
    return timeString(timeLocal, 2);   
}

/**
 * @param {Number} localtime => 通过 getLocalTime 获得的时间
 * 计算太阳高度角
 */
function calcelevation(t, localtime, latitude, longitude, zone) {
    var eqTime = calcEquationOfTime(t); // 时差
    var theta = calcSunDeclination(t); // 赤纬
    var solarTimeFix = eqTime + 4.0 * longitude - 60.0 * zone;

    var trueSolarTime = localtime + solarTimeFix;
    
    trueSolarTime = correctRange(trueSolarTime, 1440);

    var hourAngle = trueSolarTime / 4.0 - 180.0;
    if (hourAngle < -180) hourAngle += 360.0;

    var haRad = degToRad(hourAngle);
    var csz = Math.sin(degToRad(latitude)) * Math.sin(degToRad(theta)) + Math.cos(degToRad(latitude)) * Math.cos(degToRad(theta)) * Math.cos(haRad)
    if (csz > 1.0) {
        csz = 1.0
    } else if (csz < -1.0) {
        csz = -1.0
    }
    var zenith = radToDeg(Math.acos(csz))

    var exoatmElevation = 90.0 - zenith

    if (exoatmElevation > 85.0) {
        var refractionCorrection = 0.0;
    } else {
        var te = Math.tan(degToRad(exoatmElevation));
        if (exoatmElevation > 5.0) {
            var refractionCorrection = 58.1 / te - 0.07 / (te * te * te) + 0.000086 / (te * te * te * te * te);
        } else if (exoatmElevation > -0.575) {
            var refractionCorrection = 1735.0 + exoatmElevation * (-518.2 + exoatmElevation * (103.4 + exoatmElevation * (-12.79 + exoatmElevation * 0.711)));
        } else {
            var refractionCorrection = -20.774 / te;
        }
        refractionCorrection = refractionCorrection / 3600.0;
    }

    var solarZen = zenith - refractionCorrection;

    if (solarZen > 108.0) {
        document.getElementById("azbox").value = "dark";
        document.getElementById("elbox").value = "dark";
    } else {
        // 太阳高度角
        return Math.floor((90.0 - solarZen) * 100 + 0.5) / 100.0
    }
}

/**
 * @param {Number} ZeroAzimuth 零方位角
 * 零方位角 = 北, ZeroAzimuth=0
 * 零方位角 = 南, ZeroAzimuth=180
 */
function calcazimuth(t, localtime, latitude, longitude, zone, ZeroAzimuth) {
    var eqTime = calcEquationOfTime(t); // 时差
    var theta = calcSunDeclination(t); // 赤纬
    var solarTimeFix = eqTime + 4.0 * longitude - 60.0 * zone;

    var trueSolarTime = localtime + solarTimeFix;
    
    trueSolarTime = correctRange(trueSolarTime, 1440);

    var hourAngle = trueSolarTime / 4.0 - 180.0;
    if (hourAngle < -180) hourAngle += 360.0;

    var haRad = degToRad(hourAngle);
    var csz = Math.sin(degToRad(latitude)) * Math.sin(degToRad(theta)) + Math.cos(degToRad(latitude)) * Math.cos(degToRad(theta)) * Math.cos(haRad)
    if (csz > 1.0) {
        csz = 1.0
    } else if (csz < -1.0) {
        csz = -1.0
    }
    var zenith = radToDeg(Math.acos(csz))
    var azDenom = (Math.cos(degToRad(latitude)) * Math.sin(degToRad(zenith)))
    
    if (Math.abs(azDenom) > 0.001) {
        azRad = ((Math.sin(degToRad(latitude)) * Math.cos(degToRad(zenith))) - Math.sin(degToRad(theta))) / azDenom
        if (Math.abs(azRad) > 1.0) {
            if (azRad < 0) {
                azRad = -1.0
            } else {
                azRad = 1.0
            }
        }
        var azimuth = 180.0 - radToDeg(Math.acos(azRad))
        if (hourAngle > 0.0) {
            azimuth = -azimuth
        }
    } else {
        if (latitude > 0.0) {
            azimuth = 180.0
        } else {
            azimuth = 0.0
        }
    }
    if (azimuth < 0.0) {
        azimuth += 360.0
    }
    return (Math.floor(azimuth * 100 + 0.5) - ZeroAzimuth * 100) / 100.0;
}

相关方法

• 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)