A 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
}

 

[Subdot]

Calculations below to arrive at the Naval Observatory formula, if anyone is interested, and so I don't forget either! It was also fun to work through.

We start with

$$r = a\frac{1 - e^2}{1 + e\cos v},$$

where v is the true anomaly, a = 1.000001018 AU, and e is the orbital eccentricity = 0.01671.

We then use the series expansion to write the true anomaly in terms of the mean anomaly M,

$$ v = M + (2e - \frac{1}{4}e^3)\sin M + \frac{5}{4}e^2\sin 2M + \frac{13}{12}e^3\sin 3M +\mathcal{O}(e^4).$$

Then expand the earth-sun distance equation about e to get

$$r/a = 1 - e^2 - e \cos v + e^2 \cos^2 v + e^3 \cos v - e^3 \cos^3 v - e^4 \cos^2 v + e^4 \cos^4 v + \mathcal{O}(e^5),$$

but we really just need to keep terms of $\mathcal{O}(e^2)$ to reproduce the Naval Observatory's approximation. We could drop higher order terms right now, but I will leave them for the time being to see a fuller expansion. Now plug in the equation for $v(M)$ into the series expansion. The result is

$$\begin{align}r/a &= 1 - e^2 - e \cos (M + (2e - \frac{1}{4}e^3)\sin M + \frac{5}{4}e^2\sin 2M + \frac{13}{12}e^3\sin 3M) \nonumber\\
&+ e^2 \cos^2 (M + (2e - \frac{1}{4}e^3)\sin M + \frac{5}{4}e^2\sin 2M + \frac{13}{12}e^3\sin 3M) \nonumber\\
&+ e^3 \cos (M + (2e - \frac{1}{4}e^3)\sin M + \frac{5}{4}e^2\sin 2M + \frac{13}{12}e^3\sin 3M) \nonumber\\
&- e^3 \cos^3 (M + (2e - \frac{1}{4}e^3)\sin M + \frac{5}{4}e^2\sin 2M + \frac{13}{12}e^3\sin 3M)\nonumber \\
&- e^4 \cos^2 (M + (2e - \frac{1}{4}e^3)\sin M + \frac{5}{4}e^2\sin 2M + \frac{13}{12}e^3\sin 3M) \nonumber\\
&+ e^4 \cos^4 (M + (2e - \frac{1}{4}e^3)\sin M + \frac{5}{4}e^2\sin 2M + \frac{13}{12}e^3\sin 3M).\nonumber\end{align}$$

Each of the cosines can be expanded using the cosine sum identity, e.g.,

$$\begin{align}&\cos(M)\cos((2e - \frac{1}{4}e^3)\sin M + \frac{5}{4}e^2\sin 2M + \frac{13}{12}e^3\sin 3M) \nonumber\\
&- \sin(M)\sin((2e - \frac{1}{4}e^3)\sin M + \frac{5}{4}e^2\sin 2M + \frac{13}{12}e^3\sin 3M).\nonumber\\\end{align}$$

At this point, it should be apparent that to expand to $\mathcal{O}(e^2)$, the only term in the arguments for the cosine that are needed are $v = M + (2e - \frac{1}{4}e^3)\sin M$ or more shortly still $v = M + 2e\sin M$, and we only need to keep the cosines up to $\mathcal{O}(e^2)$ also. Nevertheless, proceeding with the expansion, we get the following.

$$\begin{align}r/a &= 1 - e \cos M - e^2 \cos 2 M + e^2\cos^2 M \nonumber\\
&+ (e^3 (17 \cos M - 9 \cos 3 M))/8 - 4 e^3 \cos M \sin^2 M - e^3 \cos^3 M \nonumber\\
&+ (2 e^4 (1 + 8 \cos 2 M) \sin^2 M)/3 + (e^4 (-1 - 20 \cos 2 M + 13 \cos 4 M))/8 \nonumber\\
&+ e^4(\cos^4 M + 6 \cos^2 M \sin^2 M) \nonumber\\
&+ \mathcal{O}(e^5), \nonumber\\
\nonumber\end{align}$$

To $\mathcal{O}(e^2)$, we simply just get (plugging in the values of a and e at the end),

$$\begin{align}r/a &= 1 - e \cos M - e^2 \cos 2 M + e^2\cos^2 M \nonumber\\
&= 1 - e \cos M - e^2 \cos 2 M + e^2\frac{1}{2}(1 + \cos 2M) \nonumber\\
&= 1 + \frac{e^2}{2} - e \cos M - \frac{e^2}{2} \cos 2M \nonumber\\
&\rightarrow r = 1.00014 – 0.01671 \cos M – 0.000139612 \cos 2M, \nonumber\end{align}$$

which is the result we were looking for!

As for accuracy, we estimate an upper bound on the error from the maximum value of the next higher order term, i.e., $a((e^3 (17 \cos M - 9 \cos 3 M))/8 - 4 e^3 \cos M \sin^2 M - e^3 \cos^3 M)$. The expression simplifies down to $\frac{3}{2} c^3 (\cos M - \cos^3 M)$. Its maximum value is $2.7*10^{-6}$ AU.

As for the time period when the equation is accurate, I haven't calculated (maybe I'll do it at some point; maybe not). It likely depends on when the value $e = 0.01671$ is accurate for the orbital eccentricity.

 

[D H]

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

Solar noon can happen before noon according to the clock. If they didn't subtract that half day they might well end up with the time of solar noon on the next day.

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?

You already got this in your second post, so I not answering it here.

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.

The goal is to provide the next sunrise after local midnight (approximately), the next sunset after the next sunrise (once again, approximately), and the next solar noon after local midnight. It's up to the user to know that 1 PM (local) most likely is after sunrise and solar noon.

 

[Subdot]

Thanks for the response! Sorry for the huge delay. I had a lot to do and had some travel.

Solar noon can happen before noon according to the clock. If they didn't subtract that half day they might well end up with the time of solar noon on the next day.

I'm struggling to see how this might happen. Could you walk me through it a bit? Suppose on the first pass that 1/20/2025 11:35 am local time is what is computed for solar noon. This is converted to a Julian day and then placed into the algorithm for the second pass. The time is still 11:35 am (however many minutes after local midnight that is) then when passed into the algorithm for the second pass. I don't see how it's possible to end up with solar noon on the next day, since the algorithm's input is 11:35 am 1/20/2025 local time and all quantities will be computed for that time, which is so close to local solar noon on 1/20/2025 that I would think the local solar noon for 1/20/2025 would still be picked out.

On the other hand, suppose we subtract half a day before the second pass. The input time for the second pass will now be 1/19/2025 11:25pm. This is after local midnight on 1/19/2025 but before local midnight on 1/20/2025. So wouldn't the local solar noon for 1/19/2025 then be computed?

The goal is to provide the next sunrise after local midnight (approximately), the next sunset after the next sunrise (once again, approximately), and the next solar noon after local midnight. It's up to the user to know that 1 PM (local) most likely is after sunrise and solar noon.

So if I'm understanding you here, the user has to know that their input will result in the solar times of interest, despite the timezone conversion from their local midnight to UTC. This means that NOAA is incorrect in its calculation by taking UTC midnight while still using the user's local date (i.e., NOAA assumes that midnight 1/6/2025 UTC is the correct input for the algorithm's first pass in the example of my OP), and I am correct that the local midnight for the user's date needs to be converted to UTC first?

 

[D H]

I previously gave an incorrect explanation for the subtraction of half of a day in the second pass at estimating solar noon. Here is a correct explanation.

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?

The line prior to the calculation of newt (the calculation of solNooNOffset) can be merged with the line you are questioning as

var newt = calcTimeJulianCent(jd - longitude/360.0 - eqTime/1440.0)

The -0.5 vanished in the reformulated version of the calculation of newt because the leading value in the calculation of solNoonOffset is 720 minutes; i.e., half of a day. Compare the reformulated calculation of newt with the calculation of tnoon a few lines above:

var tnoon = calcTimeJulianCent(jd - longitude/360.0)

The only difference between the calculations of newt and tnoon a few lines above are the subtraction of eqTime/1440.0 in the calculation of newt.

I am more concerned with NOAA naming that variable tnoon compared to their subtracting 0.5 and then adding 720/1440 (=0.5). The time jd - longitude/360.0 is local solar mean midnight expressed in UTC.

 

[Subdot]

I previously gave an incorrect explanation for the subtraction of half of a day in the second pass at estimating solar noon. Here is a correct explanation.

Thank you! That is indeed what is happening. So now the big question: Is this correct? They are calculating the equation of time at the estimated local midnight in UTC for both the first and second pass. This is deliberate: if they wanted to calculate the equation of time at the estimated solar noon, they would not have subtracted the half day.

I would have thought the goal was to calculate the equation of time as close to the sun's actual position at the solar time of interest. I'm therefore inclined to think this subtraction of half a day is incorrect, and we would want to plug into the second pass the result of the first pass--our initial estimate of solar noon (just like is done for sunrise/sunset). Am I missing something?


I wonder if this is related to my question about which julian day to pick. I'll have to think on it more, unless you or someone else figures it out.

 

[WernerQH]

They are calculating the equation of time at the estimated local midnight in UTC for both the first and second pass. This is deliberate: if they wanted to calculate the equation of time at the estimated solar noon, they would not have subtracted the half day.

I don't think so. They are strictly adhering to the definition of the julian date:

The number of days since noon on January 1, 4713 BC, plus the decimal fraction of a day that has elapsed since the preceding noon.

The Julian date (JD) of any instant is the Julian day number plus the fraction of a day since the preceding noon in Universal Time. Julian dates are expressed as a Julian day number with a decimal fraction added.

In the old days it was considered helpful (at least for European astronomers) not to have a change of the day number during the nightly observations.

 

[Subdot]

I don't think so. They are strictly adhering to the definition of the julian date:


In the old days it was considered helpful (at least for European astronomers) not to have a change of the day number during the nightly observations.

Good thought, but I don't see how this is what is going on. The user inputs 1/6/2025 (any timezone). The julian day is computed for 1/6/2025 with the function getJD() in the OP. 1524.5 is subtracted at the end of this conversion, so the julian day is for 1/6/2025 at UTC midnight (so I understand). This is the jd that is the input for the sunrise, sunset, and solar noon functions.

The input jd must then be at midnight (which is then converted to local mean midnight by subtracting the longitude), not noon, so the reasoning of my previous posts holds that the second pass is calculated for local mean midnight.

Did I miss something?

 

[WernerQH]

The user inputs 1/6/2025 (any timezone). The julian day is computed for 1/6/2025 with the function getJD() in the OP. 1524.5 is subtracted at the end of this conversion, so the julian day is for 1/6/2025 at UTC midnight (so I understand).

The point of subtracting 1524.5 is to conform to the conventional definition. It produces the julian date for January 6, 2025, 0:00. (Or perhaps January 7? I'm not sure.)

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?

And they subtract another half a day to get the time for noon. Maybe they should add half a day? There is a "modified julian date" (MJD) that is the julian date with 2400000.5 subtracted. If I have looked it up correctly, it should be MJD = 60681 for Jan 6, 2025.

 

[WernerQH]

2025-01-06 0:00 UT = MJD 60681 = JD 2 460 681.5
This is what getJD(2025,1,6) returns. It seems noon, sunrise and sunset are computed for the previous day.

 

[Subdot]

2025-01-06 0:00 UT = MJD 60681 = JD 2 460 681.5
This is what getJD(2025,1,6) returns. It seems noon, sunrise and sunset are computed for the previous day.

Confirming after just playing around with codepen that after choosing 1/6/2025, 2460681.5 is the input for solar noon, sunrise, and sunset, which means the time is computed at midnight, not noon.

By various calculators, including the JPL, converting the 2460681.5 to a calendar date is not the previous day though: it is still 1/6/2025 at midnight: https://ssd.jpl.nasa.gov/tools/jdc/#/cd

Edit 2: Unless you were referring to the julian day? In that case, noon for 1/6/2025 after midnight on 1/6/2025 would be 2460682. Does the algorithm then look for solar noon for 1/5/2025 because the Julian day input is still within the 2460681 day? If so, then yeah, I guess it is calculating noon for the previous day. However, if the algorithm looks for solar noon after 2460681.5, then the next noon is 2460682, so it would not be the previous day but 1/6/2025 noon as desired. Based on DH's earlier comments in the thread, I would think that the latter is the behavior: solar noon after local mean midnight is found.

I couldn't find an independent calculator to check it, but the calcJDFromJulianCent() function in the NOAA calclator converts the julian century (0.2501493649707162 on the first pass) back to the original julian day (2460681.7055555554 = 2460681.5 - longitude info). So there is likely no rounding that happens during the julian century conversion that would change the input to a time at noon.

The javascript function to convert to julian century from julian day would indicate no rounding either: it is double-arithmetic.

function calcTimeJulianCent(jd) {
    var T = (jd - 2451545.0)/36525.0
    return T
}

function calcJDFromJulianCent(t) {
    var JD = t * 36525.0 + 2451545.0
    return JD
}

Edit: There are intriguing comments on NOAA's old solar calculator, which suggest they at least originally were trying to compute sunrise, sunset, and solar noon using local mean solar noon on the first pass; not sure what to make of the second pass, which is the one that is under question in this thread, though at a glance it appears they subtracted the half day to get back to midnight so that the second pass receives an estimated solar noon, instead of midnight (I don't know if they actually succeeded in doing that if that was their goal; haven't looked at the rest of the code).

(Note that their old calculator uses the older convention of lat/long, which is opposite of what is used today and their new calculator; hence they add longitude instead of subtract).

[CODE lang="javascript" highlight="13-15, 46-48, 53, 97-99"]
//***********************************************************************/
//* Name: calcSolNoonUTC */
//* Type: Function */
//* Purpose: calculate the Universal Coordinated Time (UTC) of solar */
//* noon for the given day at the given location on earth */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* longitude : longitude of observer in degrees */
//* Return value: */
//* time in minutes from zero Z */
//***********************************************************************/

function calcSolNoonUTC(t, longitude)
{
// First pass uses approximate solar noon to calculate eqtime
var tnoon = calcTimeJulianCent(calcJDFromJulianCent(t) + longitude/360.0);
var eqTime = calcEquationOfTime(tnoon);
var solNoonUTC = 720 + (longitude * 4) - eqTime; // min

var newt = calcTimeJulianCent(calcJDFromJulianCent(t) -0.5 + solNoonUTC/1440.0);

eqTime = calcEquationOfTime(newt);
// var solarNoonDec = calcSunDeclination(newt);
solNoonUTC = 720 + (longitude * 4) - eqTime; // min

return solNoonUTC;
}

//***********************************************************************/
//* Name: calcSunriseUTC */
//* Type: Function */
//* Purpose: calculate the Universal Coordinated Time (UTC) of sunrise */
//* for the given day at the given location on earth */
//* Arguments: */
//* JD : julian day */
//* latitude : latitude of observer in degrees */
//* longitude : longitude of observer in degrees */
//* Return value: */
//* time in minutes from zero Z */
//***********************************************************************/

function calcSunriseUTC(JD, latitude, longitude)
{
var t = calcTimeJulianCent(JD);

// *** Find the time of solar noon at the location, and use
// that declination. This is better than start of the
// Julian day

var noonmin = calcSolNoonUTC(t, longitude);
var tnoon = calcTimeJulianCent (JD+noonmin/1440.0);

// *** First pass to approximate sunrise (using solar noon)

var eqTime = calcEquationOfTime(tnoon);
var solarDec = calcSunDeclination(tnoon);
var hourAngle = calcHourAngleSunrise(latitude, solarDec);

var delta = longitude - radToDeg(hourAngle);
var timeDiff = 4 * delta; // in minutes of time
var timeUTC = 720 + timeDiff - eqTime; // in minutes

// alert("eqTime = " + eqTime + "\nsolarDec = " + solarDec + "\ntimeUTC = " + timeUTC);

// *** Second pass includes fractional jday in gamma calc

var newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440.0);
eqTime = calcEquationOfTime(newt);
solarDec = calcSunDeclination(newt);
hourAngle = calcHourAngleSunrise(latitude, solarDec);
delta = longitude - radToDeg(hourAngle);
timeDiff = 4 * delta;
timeUTC = 720 + timeDiff - eqTime; // in minutes

// alert("eqTime = " + eqTime + "\nsolarDec = " + solarDec + "\ntimeUTC = " + timeUTC);

return timeUTC;
}

//***********************************************************************/
//* Name: calcSunsetUTC */
//* Type: Function */
//* Purpose: calculate the Universal Coordinated Time (UTC) of sunset */
//* for the given day at the given location on earth */
//* Arguments: */
//* JD : julian day */
//* latitude : latitude of observer in degrees */
//* longitude : longitude of observer in degrees */
//* Return value: */
//* time in minutes from zero Z */
//***********************************************************************/

function calcSunsetUTC(JD, latitude, longitude)
{
var t = calcTimeJulianCent(JD);

// *** Find the time of solar noon at the location, and use
// that declination. This is better than start of the
// Julian day

var noonmin = calcSolNoonUTC(t, longitude);
var tnoon = calcTimeJulianCent (JD+noonmin/1440.0);

// First calculates sunrise and approx length of day

var eqTime = calcEquationOfTime(tnoon);
var solarDec = calcSunDeclination(tnoon);
var hourAngle = calcHourAngleSunset(latitude, solarDec);

var delta = longitude - radToDeg(hourAngle);
var timeDiff = 4 * delta;
var timeUTC = 720 + timeDiff - eqTime;

// first pass used to include fractional day in gamma calc

var newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440.0);
eqTime = calcEquationOfTime(newt);
solarDec = calcSunDeclination(newt);
hourAngle = calcHourAngleSunset(latitude, solarDec);

delta = longitude - radToDeg(hourAngle);
timeDiff = 4 * delta;
timeUTC = 720 + timeDiff - eqTime; // in minutes

return timeUTC;
}
[/CODE]

 

[WernerQH]

Confirming after just playing around with codepen that after choosing 1/6/2025, 2460681.5 is the input for solar noon, sunrise, and sunset, which means the time is computed at midnight, not noon.

Yes, the input time is consistently midnight.

Sorry, I was confused and hadn't looked at the code closely enough. In the line preceding the one you first quoted, a variable solNoonOffset is defined

var solNoonOffset = 720.0 - (longitude * 4) - eqTime // in minutes

where an offset of half a day (720 minutes) is added, so everything seems to be okay.

The javascript function to convert to julian century from julian day would indicate no rounding either: it is double-arithmetic.

It's just a conversion of the same time from days to centuries (with a shift of origin to the standard epoch J2000.0), presumably for the calculation of the equation of time.

 

[Subdot]

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 now understand why my naive way of doing things here is incorrect and NOAA's way is correct, thanks to the discussion. Basically, my concerns about finding the Julian day for the user's local midnight are already taken care of by subtracting the longitude in the solar noon function var tnoon = calcTimeJulianCent(jd - longitude/360.0). If we convert the user's local midnight to UTC first, this subtraction will overcorrect and result in the incorrect date. However, this then leaves me mystified for the sunrise/sunset equation, which does not calculate at the user's local midnight.

Using my example above, midnight user's local time in China is 1/6/2025 12am (UTC+8), which corresponds to 1/5/2025 4pm UTC. Following NOAA's calculator, we start with 1/6/2025 12am UTC. We then subtract the longitude to get local mean solar midnight, which results in 2460681.2105833334 -> 1/5/2025 3:03 pm UTC, which is an alright approximation to the user's local midnight (results in 1/6/2025 11:03 pm (UTC+8)).

If instead, we converted the date to UTC first (1/5/2025 4pm), subtracting the longitude would change this to 1/5/2025 7am, which is 1/5/2025 3pm (UTC+8)--not anywhere close to the local midnight.

We could use an example east of UTC too and see that the date is correct. This is more believable, since timezones east of UTC will have a midnight that clearly is on the same date as the UTC date (cause you add the timezone offset, not subtract).



The reason it works is that UTC is always +/-12 < 24 hours of the different timezones of the world (more or less; we have +14 but that's still less than 24), and the timezones roughly correspond to local mean solar midnight, so the local mean solar midnight always arrives at the midnight of the correct local date that the user requests.

 

[Subdot]

It seems to me, the way to determine whether NOAA's approximation is better than alternatives I've suggested in the thread is to know:

1) Is it true that we want the input time for the second pass to the equation of time, declination, etc. to be as close to the actual time as we can get? I think this is true, but I don't know for sure. Confirmation would be helpful before I ditch NOAA for doing my own thing (and maybe email NOAA, though they probably won't respond)!

2) If the input for the first pass is not approximately local midnight or later (like local noon), is it possible for the equation of time, declination, etc. to compute for the previous day, resulting in something inaccurate? It could be that what NOAA does on the first pass for sunrise/sunset (UTC midnight, not local midnight) doesn't change anything significant. This is something I don't even have a clue for the answer. I would need to understand the equation of time, declination, and so forth better than I do. My understanding is that they are just simple time dependent equations, so the closer the time is to the previous day's sunrise/sunset, etc., the less accurate the approximation will be (because it is based on the sun's position closer to the previous day's sunrise/sunset, rather than closer to the desired day's sunrise/sunset).

Based on the old calculator and its comments, I suspect that an error was made when moving the code to the new calculator, though even the old calculator for solar noon used local midnight for its first pass, despite calling the variable tnoon still. Answering these questions would help me see whether the code is mistaken/not as accurate as it could be or maybe once was.

@WernerQH @D H

 

[WernerQH]

1) Is it true that we want the input time for the second pass to the equation of time, declination, etc. to be as close to the actual time as we can get? I think this is true, but I don't know for sure.

Yes, of course one should calculate the relevant parameters for the exact time at which the sun rises or sets at one's location. In principle, this could involve even several iterations. But it would be silly to aim for millisecond accuracy. Earth's rotation is not so regular. Today's difference between UT1 and UTC, for example, is about 48 milliseconds.

2) If the input for the first pass is not approximately local midnight or later (like local noon), is it possible for the equation of time, declination, etc. to compute for the previous day, resulting in something inaccurate? It could be that what NOAA does on the first pass for sunrise/sunset (UTC midnight, not local midnight) doesn't change anything significant.

You can express times in UT1, UTC, TAI, or local time. There are even several kinds of ephemeris time, differing by seconds, if I remember correctly. The calculations can be done using any of those, and converted to local time in the end. The result should always be the same. But of course the best initial estimate for local noon is 12 hours after local midnight, or (720 - longitude*4) minutes after midnight UTC.

The declination of the sun can change by up to 0.4 degrees per day, or 6 minutes of arc within 6 hours, so this can have a slight effect on the time of sunset. Likewise, the equation of time can vary by up to 30 seconds per day, so the time of sunset might be off by 7 seconds, if you use the equation of time calculated for local noon. But I think this is neglected in most cases.

A look at the Wikipedia page might be helpful.

 

[Subdot]

You can express times in UT1, UTC, TAI, or local time. There are even several kinds of ephemeris time, differing by seconds, if I remember correctly. The calculations can be done using any of those, and converted to local time in the end. The result should always be the same. But of course the best initial estimate for local noon is 12 hours after local midnight, or (720 - longitude*4) minutes after midnight UTC.

The declination of the sun can change by up to 0.4 degrees per day, or 6 minutes of arc within 6 hours, so this can have a slight effect on the time of sunset. Likewise, the equation of time can vary by up to 30 seconds per day, so the time of sunset might be off by 7 seconds, if you use the equation of time calculated for local noon. But I think this is neglected in most cases.

A look at the Wikipedia page might be helpful.

Thank you! Interesting that the Wikipedia page subtracts/adds the hour angle from the time of true solar transit, rather than from (720 - longitude*4). I would naively think that the Wikipedia page's suggestion is in theory more accurate? I tried an example, and this produces significant differences up to several seconds (mathematically, cause the eqTime is calculated at noon, rather than at sunrise/sunset, so off by 6 hour -> the 7 seconds you mentioned).

Yet NOAA can't be that far off: as I recall, it matches timeanddate.com to 3 seconds or so for sunrise/sunset (NOAA is more different for solar noon, but we know why now); unless timeanddate.com is also just far off. The case I just tried shows that using Wikipedia's suggestion makes the calculation different in the wrong direction (1/6/2025 at 40, -74 from timeanddate.com says 7:17 am for sunrise, while I am getting 7:18:06 am using Wikipedia's suggestion and 7:18:00 am using NOAA's current way of doing things).

Mathematically, Wikipedia's suggestion (expressed using NOAA's definitions) is to use

sunrise/sunset = 720 - 4*longitude - 4*hourangle - eqTime(evaluated at approximately true solar noon),

whereas NOAA is calculating sunrise/sunset as

sunrise/sunset = 720 - 4*longitude - 4*hourangle - eqTime(evaluated at approximately true sunrise/sunset)


Who is correct or more correct? I guess it depends on whether the hour angle is defined relative to the mean solar times or to the true solar times. Wikipedia's definition of the hour angle suggests it is defined relative to true solar times, which leaves me mystified for the reasons I noted above.


On the first pass: it doesn't matter to seconds accuracy

From playing around a bit with the algorithm, thanks to the second pass, it turns out that it doesn't matter if the first pass is midnight at the UTC timezone (which is what I intended when I said UTC midnight in my previous post; I did not mean midnight expressed in UTC; apologies for not speaking more clearly), i.e., jd, or if it is local midnight (jd - longitude*4), or if it is local noon (jd - longitude*4 + 720): any of these initial guesses results in the same value to milliseconds (tenth of a millisecond for sunrise, for the date and location that I played with) after the second pass.

I would assume that the best initial guess, based on our discussion, is local mean solar noon for all three times (thinking of approximate sunrise as noon - 6 hours and approximate sunset as noon + 6 hours), just like NOAA did in their old calculator. Local mean solar midnight should on average be the same accuracy for sunrise as for local mean solar noon, though it was different (to tenths of a millisecond, if I recall correctly) in the cases I examined. So I'll go with that for conceptual reasons, even if it doesn't affect the end result due to the second pass.

 

[WernerQH]

Mathematically, Wikipedia's suggestion (expressed using NOAA's definitions) is to use

sunrise/sunset = 720 - 4*longitude - 4*hourangle - eqTime(evaluated at approximately true solar noon),

whereas NOAA is calculating sunrise/sunset as

sunrise/sunset = 720 - 4*longitude - 4*hourangle - eqTime(evaluated at approximately true sunrise/sunset)


Who is correct or more correct? I guess it depends on whether the hour angle is defined relative to the mean solar times or to the true solar times. Wikipedia's definition of the hour angle suggests it is defined relative to true solar times, which leaves me mystified for the reasons I noted above.

The Wikipedia algorithm is simplified -- the hour angles for sunrise and sunset are assumed equal, whereas during winter and spring, when the sun's declination increases and the days are getting longer, the afternoons are slightly longer than the mornings, and vice versa in summer and autumn. The hour angle is a local, purely geometric quantity: the angle between the meridian (south) and the point of sunrise/sunset, as seen from the celestial north pole. So the NOAA algorithm is more accurate. What I found a little confusing is that the sunrise/sunset times (which of course depend on location) are returned as UTC values. But it makes sense because equation of time and declination are recalculated, and these functions do not depend on location.