Understanding NOAA and Naval Observatory's Solar Calculations

[Subdot]

TL;DR Summary

Which approximations for calculating sunrise, sunset, solar noon are most accurate? NOAA and Naval Observatory do something different for earth-sun distance: where does this equation come from? NOAA also subtracts half a day in its algorithm for solar noon but not sunrise/set: why? Also, why convert to Julian day before converting to UTC?

I have three questions in relation to algorithms used in calculating sunrise, sunset, and solar noon.

I'm looking at the NOAA javascript (view source on their solar calculator page and then find the main.js file) and comparing some of the intermediate calculations to other algorithms, such as in the Naval Observatory's earth-sun distance approximation. I have noticed some differences, and I want to know where these approximations are coming from and what is the error in using the different methods? I do not have access to the Astronomical Almanac or other paywalled references. I saw some discussion on some of these questions in an earlier physicsforums thread, the answer to some of my questions appears to be that several equations are being smooshed together and expanded, though I don't know which ones.

I also have a basic question about the use of the Julian day in the calculations for which I am having difficulty finding an answer online.


1) Why does the NOAA subtract a half day in the second pass at estimating solar noon?

Other calculators just do like NOAA did for sunrise: add the time for the first guess to the date at midnight; NOAA does this but also subtracts half a day.

jd below is julian day; timezone is the timezone offset to UTC. The rest is self-explanatory.

I'm looking at the line that goes:

var newt = calcTimeJulianCent(jd - 0.5 + solNoonOffset/1440.0)

See how they subtract half a day before using newt as input for the equation of time?

Here is the full function for solar noon.

function calcSolNoon(jd, longitude, timezone) {
    var tnoon = calcTimeJulianCent(jd - longitude/360.0)
    var eqTime = calcEquationOfTime(tnoon)
    var solNoonOffset = 720.0 - (longitude * 4) - eqTime // in minutes
    var newt = calcTimeJulianCent(jd - 0.5 + solNoonOffset/1440.0)
    eqTime = calcEquationOfTime(newt)
    var solNoonLocal = 720 - (longitude * 4) - eqTime + (timezone*60.0)// in minutes
    while (solNoonLocal < 0.0) {
        solNoonLocal += 1440.0;
    }
    while (solNoonLocal >= 1440.0) {
        solNoonLocal -= 1440.0;
    }

    return solNoonLocal
}

...And here is sunrise/sunset.

// rise = 1 for sunrise, 0 for sunset
function calcSunriseSet(rise, JD, latitude, longitude, timezone) {

    var timeUTC = calcSunriseSetUTC(rise, JD, latitude, longitude);
    var newTimeUTC = calcSunriseSetUTC(rise, JD + timeUTC/1440.0, latitude, longitude);
    if (isNumber(newTimeUTC)) {
        var timeLocal = newTimeUTC + (timezone * 60.0)
        var riseT = calcTimeJulianCent(JD + newTimeUTC/1440.0)
        var riseAzEl = calcAzEl(riseT, timeLocal, latitude, longitude, timezone)
        var azimuth = riseAzEl.azimuth
        var jday = JD
        if ( (timeLocal < 0.0) || (timeLocal >= 1440.0) ) {
            var increment = ((timeLocal < 0) ? 1 : -1)
            while ((timeLocal < 0.0)||(timeLocal >= 1440.0)) {
                timeLocal += increment * 1440.0
                jday -= increment
            }
        }

    } else { // no sunrise/set found

        var azimuth = -1.0
        var timeLocal = 0.0
        var doy = calcDoyFromJD(JD)
        if ( ((latitude > 66.4) && (doy > 79) && (doy < 267)) ||
             ((latitude < -66.4) && ((doy < 83) || (doy > 263))) ) {
            //previous sunrise/next sunset
            jday = calcJDofNextPrevRiseSet(!rise, rise, JD, latitude, longitude, timezone)
        } else {   //previous sunset/next sunrise
            jday = calcJDofNextPrevRiseSet(rise, rise, JD, latitude, longitude, timezone)
        }
    }

    return {"jday": jday, "timelocal": timeLocal, "azimuth": azimuth}
}

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;    // in minutes

    return timeUTC
}

2) Where does the Naval observatory's equation for approximating the earth-sun distance come from (UPDATE: I have now been able to reproduce their results)? Is it more accurate than using the polar ellipse equation?

The Naval observatory chooses to use the following equation to approximate the earth-sun distance.

$$R = 1.00014 – 0.01671 \cos g – 0.00014 \cos 2g,$$

where g is the mean anomaly of the Sun and R is the earth-sun distance.

Where does this equation come from? Is it coming from the Jet Propulsion lab, e.g., DE440, and the chebyshev polynomials (but they seem based on Julian days, not the anomaly directly)? It does not appear to be a series expansion of the equation the NOAA uses to calculate the earth-sun distance. Basically, the NOAA is using $r = a(1-e^2)/(1 + e\cos(v))$, where v is the true anomaly and e is the earth's orbital eccentricity, so it could be there is another expansion involved here to write NOAA's expression in terms of the mean anomaly. Based on the earlier physicsforums thread, it seems my suspicion is correct that the mean anomaly is being plugged into the equation for the true anomaly, though I haven't been able to reproduce the Naval Observatory's formula yet.

UPDATE: The equation is indeed coming from plugging in the equation for the true anomaly in terms of the mean anomaly and then expanding about the eccentricity parameter, taking it to be a constant $e = 0.01671$. My question about accuracy remains then because...why? Why make this approximation when the orbital trajectory could be calculated without doing so many series expansions and when the orbital eccentricity could be found as a function of time? I suspect NOAA's way may be more accurate, but I am seeking confirmation!

My other question on this though is: What is the error in the formula from the Naval observatory? What is the error in the formula for the NOAA? Which one is more accurate? If the Naval observatory is just some expansion of NOAA's equation, it would seem that NOAA's equation is more accurate?

Here is NOAA's function for calculating the earth-sun distance for reference (which as I recall is also what Meeus says to do), where t is in julian centuries.

function calcSunRadVector(t) {
    var v = calcSunTrueAnomaly(t);
    var e = calcEccentricityEarthOrbit(t);
    var R = (1.000001018 * (1 - e * e)) / (1 + e * Math.cos(degToRad(v)));
    return R;        // in AUs
}

3) Why does NOAA not convert the user's input date to UTC before converting to the Julian day?

If the user inputs 1/6/2025 1pm (UTC+8), the Julian day is computed for 1/6/2025 with no time component (so basically, computing the Julian day for 1/6/2025 at midnight UTC). This is the Julian day that serves as input to the sunrise/set and solar noon functions.

However, I would naively think the user's input date at midnight should be converted to UTC first 1/6/2025 at midnight (UTC+8) --> 1/5/2025 4pm UTC, and so the julian day should be computed for 1/5/2025 4pm.

I've run through some scenarios and can see that somehow the NOAA way of doing things works, but I don't understand why. What is going on conceptually here? Why is this better? Why would converting the user's input date at midnight to UTC first as I would naively think be incorrect or not as good an approximation?


For reference, here is the function NOAA uses to convert a date to Julian day. Pretty standard stuff though.

function getJD(year, month, day) {
    if (month <= 2) {
        year -= 1
        month += 12
    }
    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
}