solar – Sunrise and Sunset¶
This module computes sunrise and sunset.
See https://en.wikipedia.org/wiki/Sunrise_equation
Finally, see https://gml.noaa.gov/grad/solcalc/sunrise.html which seems to produce different, more accurate results. The associated Excel spreadsheet seems more useful because it seems to be the preferred approach and provides test data.
Supplemental information: https://gml.noaa.gov/grad/solcalc/calcdetails.html
Notation:
- \(\phi\)
Latitude, north is positive.
- \(\lambda\)
Longitude, east is positive.
- \(t\)
Time zone offset in hours. Converts UTC to local time, and accounts for savings vs. daylight time.
The procedure is as follows:
D is the given date as an ordinal number after 1900. The spreadsheet representation of dates provides a helpful approximation.
date - (1900-1-1) + 2.E is the time past midnight of D as a float, in minutes.
F is the Julian day number. (\(t\) is the timezone offset, in hours.)
\[F = D + 2415018.5 + E - \frac{t}{24}\]G is the Julian century.
\[G = \frac{F - 2451545}{36525}\]I is the sun’s geometric mean longitude (in degrees).
\[I = (280.46646 + 36000.76983 G + 0.0003032 G^2) \mod 360\]J is the sun’s geometric mean anomaly (in degrees).
\[J = 357.52911 + 35999.05029 G - 0.0001537 G^2\]K is the eccentricity of earth’s orbit.
\[K = 0.016708634 - 0.000042037 G + 0.0000001267 G^2\]L is the sun’s equation of center.
\[L = (1.914602 - 0.004817 G + 0.000014 G^2) \sin \frac{J \pi}{180} + (0.019993 - 0.000101 G) \sin \frac{2 J \pi}{180} + 0.000289 \sin \frac{2 J \pi}{180}\]M is the sun’s true longitude (in degrees), \(M = I + L\).
N is the sub’s true anomaly (in degrees), \(M = J + L\).
(Not required, but is in the original spreadsheet.) O is the sun’s radius vector (in Astronomical Units).
\[O = \frac{1.000001018 (1 - K^2)}{1 + K \cos N}\]P is the sun’s apparent longitude (in degrees).
\[P = M - 0.00569 - 0.00478 \sin \left(\frac{(125.04 - 1934.136 G)\pi}{180}\right)\]Q is the mean obliquity to the ecliptic (in degrees).
\[Q = 23 + \frac{26}{60} + \frac{21.448 - 46.815 G + 0.00059 G^2 - 0.001813 G^3}{3600}\]R is the obliquity correction (degrees).
\[R = Q + 0.00256 \cos \left(\frac{(125.04 - 1934.136 G) \pi}{180} \right)\]S is the sun’s right ascension (degrees).
\[S = \frac{180}{\pi} \arctan \left( \frac{\cos R \sin P} {\cos P} \right)\]T is the sun’s declination (degrees).
\[T = \frac{180}{\pi} \arcsin \left(\sin R \sin P \right)\]U is the “variable Y”. This is referenced as \(y\) in several variations. See https://ui.adsabs.harvard.edu/abs/1989MNRAS.238.1529H/abstract for an example.
\[U = \tan^2 \frac{\frac{R}{2}\pi}{180} = \tan^2 \frac{\pi R}{360}\]V is the “Equation” of time (in minutes) how apparent time equates to measured time. (Note, \(I\), and \(J\) are in degrees, and the radians conversions are omitted.)
\[V = 4 U \sin 2I - 2K\sin J + 4KU\sin J \cos 2I - 0.5 U^2 \sin 4I - 1.25 K^2 \sin 2J\]W is the Local Hour Angle of Sunrise (degrees). The baseline for a visible sun is \(90^{\circ}50^{\prime} = 90^{\circ}\!.833\). Other values account for civil, nautical, or astronomial twilight. (See the Twilight Details.)
\[W = \arccos \left(\frac{\cos \frac{90^{\circ}\!.833\pi}{180}}{\cos \phi \cos T} - \tan \phi \tan T\right)\]X is Solar Noon in Local Standard Time. Note that \(1440 = 24 \times 60\), the number of minutes in a day. The equation of time, \(V\), is in minutes.
\[X = \frac{720 - 4 \lambda - V + \frac{t}{60}}{1440}\]Y is sunrise.
\[Y = \frac{X - 4 W}{1440}\]Z is sunset.
\[Z = \frac{X + 4 W}{1440}\]
Twilight Details¶
Computing nautical twilight is Solar zenith angle is 102°, solar elevation angle is -12°. This is 12° beyond the horizon.
For some information, see http://www.stargazing.net/kepler/sunrise.html
This appears to be an offset to the Hour Angle of sunrise. The following uses 90.833 as the horizon with a correction for refraction of the atmosphere:
This means Nautical twilight is:
Alternative 1¶
See https://gml.noaa.gov/grad/solcalc/solareqns.PDF for a simpler (but less accurate) approach.
The above solareqns.PDF uses what appear to be deprecated Spencer equations https://www.mail-archive.com/sundial@uni-koeln.de/msg01050.html. It appears these old Fourier Transform approximations are no longer considered accurate enough. Also, this paper seems to have a number of small errors in it. (See “cost”, for example, which should be “cos”.)
Inputs:
- \(\phi\)
latitude of observer (north is positive; south is negative)
- \(\lambda\)
longitude of observer (east is positive; west is negative)
- \(N\)
Days after start of year.
- h, m, s
Current local time.
- tz
Timezone Hours offset from UTC. (For example, US_MST = -7).
Compute fraction of year, \(\gamma\), from year, \(y\), ordinal day, \(N\), and hour, \(h\).
\[\begin{split}\begin{align} d(y)& = \begin{cases} 366& \textbf{if $y \mod 4 = 0 \land (y \mod 100 \neq 0 \lor y \mod 400 = 0)$}\\ 365& \textbf{otherwise} \end{cases}\\ \gamma& = \frac{2 \pi}{d(y)}\left( N - 1 + \frac{h-12}{24}\right) \end{align}\end{split}\]
From \(\gamma\), estimate the equation of time, \(EqT\) (in minutes).
\[\begin{split}\begin{align} EqT& = 229.18 (0.000075 \\ &+ 0.001868 \cos \gamma - 0.032077 \sin \gamma\\ &- 0.014615 \cos 2\gamma - 0.040849 \sin 2\gamma) \end{align}\end{split}\]From \(\gamma\), estimate the solar declension angle, \(\delta\) (in radians).
\[\begin{split}\begin{align} \delta& = 0.006918 - 0.399912 \cos \gamma + 0.070257 \sin \gamma \\ &- 0.006758 \cos 2 \gamma + 0.000907 \sin 2 \gamma \\ &- 0.002697 \cos 3 \gamma + 0.00148 \sin 3 \gamma \end{align}\end{split}\]Find the time offset, in minutes, given the longitude, \(\lambda\), and the timezone hours, \(tz\), from UTC. (example US_MST = -7)
\[t = EqT + 4 \lambda - 60 tz\]Find the true solar time, in minutes.
\[T = 60 h + m + \frac{s}{60} + t\]Compute the Hour Angle at the Zenith. For the special case of sunrise or sunset, the zenith is set to \(90^{\circ}\!.833\) (the approximate correction for atmospheric refraction at sunrise and sunset, and the size of the solar disk), and the hour angle becomes:
\[\cos H = \frac{\cos 90.833}{\cos \phi \cos \delta}-(\tan \phi \tan \delta)\]Sunrise, \(T_r\).
\[720-4(\lambda + H) - EqT\]Sunset, \(T_s\).
\[720-4(\lambda - H) - EqT\]Noon, \(T_n\).
\[720-4\lambda - EqT\]
Alternative 2¶
See https://edwilliams.org/sunrise_sunset_algorithm.htm
This is a reference to the Almanac for computers, 1990 edition.
The official citation:
United States Naval Observatory. Nautical Almanac Office. (19801991). Almanac for computers. Washington, D.C.: Nautical Almanac Office, United States Naval Observatory.
This has a section, Sunrise, Sunset and Twilight, which we reproduce much of.
Day of the Year¶
Preliminary Information from Page B1 and B2.
- \(N\)
Day of year, integer; the time in days since the beginning of the year Range is 1 to 365 for non-leap years, 1 to 366 for leap years.
Where \(N\) is the day of the year, \(K\) is the year, \(M\) is the month (\(1 \leq M \leq 12\)), and \(I\) is the day of the month (\(1 \leq I \leq 31\)).
This is valid for any year except centurial years not evenly divisible by 400. This is valid for 2000, but not for 1900 or 2100.
Sunrise, Sunset and Twilight¶
This starts on Page B5.
For locations between latitudes 65° North and 65° South, the following algorithm provides times of sunrise, sunset and twilight to an accuracy of \(\pm 2 ^m\), for any date in the latter half of the twentieth century.
Notation:
- \(\phi\)
latitude of observer (north is positive; south is negative)
- \(\lambda\)
longitude of observer (east is positive; west is negative)
- \(M\)
Sun’s mean anomaly
- \(L\)
Sun’s true longitude
- \(RA\)
Sun’s right ascension
- \(\delta\)
Sun’s declination
- \(H\)
Sun’s local hour angle
- \(z\)
Sun’s zenith distance for sunrise or sunset. One of the following:
Zenith Choices¶ Phenomenon
z
cos z
Sunrise Sunset
90°50′
-0.01454
Civil Twilight
96°
-0.10453
Nautical Twilight
102°
-0.20791
Astronomical Twilight
108°
-0.30902
Noon
0°
+1.00000
- \(t\)
approximate time of phenomenon in days since 0 January, \(O^h\) UT
- \(T\)
local mean time of phenomenon
- \(T_U\)
universal time of phenomenon
Formulas:
Procedure:
With an initial valueof \(t\), compute \(M\) from eq. (1) and then \(L\) from eq. (2). If a morning phenomenon (sunrise or the beginning of morning twilight) is being computed, construct an initial value of \(t\) from the formula
\[t=N+(6^{h}-\lambda)/24\]Where \(N\) is the day of the year (see calendar formulas on page B1) and \(\lambda\) is the observer’s longitude expressed in hours.
If an evening phenomenon is being computed, use
\[t=N+(18^{h}-\lambda)/24\]For transit of the local meridian (i.e., noon), use
\[t=N+(12^{h}-\lambda)/24\]Solve eq. (3) for \(RA\), nothing that \(RA\) is in the same quadrant as \(L\). Transform \(RA\) to hours for later use in eq. (6).
Solve eq. (4) for \(\sin \delta\), which appears in eq. (5); \(\cos \delta\), which also is required in eq. (5), should be determined from \(\sin \delta\). While \(\sin \delta\) may be positive or negative, \(\cos \delta\) is always positive.
Solve eq. (5) for \(H\). Since computers and calculators normally give arccosine in the range 0°-180°, the correct quadrant for \(H\) can be selected according to the following rules:
rising phenomena: \(H = 360^{\circ} - \arccos x\);
setting phenomena: \(H = \arccos x\).
In other words, for rising phenomena, \(H\) must be in quadrant 3 or 4 (depending on the sign of \(\cos H\)), whereas \(H\) must be either in quadrant 1 or 2 for setting phenomena. Convert \(H\) from degrees to hours for use in eq. (6).
Compute \(T\) from eq. (6), recalling that \(H\) and \(RA\) must be expressed in hours. If \(T\) is negative or greater than \(24^h\), it should be converted to the range \(0^h - 24^h\) by adding or subtracting multiples of \(24^h\).
Compute \(T_U\) from eq. (7), where \(\lambda\) must be expression in hours. \(T_U\) is an approximation to the time of the desired rising or setting phenomenon, referred to the Greenwich meridian. If \(T_U\) is greater than \(24^h\), the phenomenon occurs on the following day, Greenwich time. If \(T_U\) is negative, the phenomenon occurs on the previous day day, Greenwich time.
Under certain conditions, eq. (5) will yield a value of \(\lvert{\cos H}\rvert > 1\), indicating the absence of the phenomenon on that day. At far northern latitudes, for example, there is continuous illumination during certain summer days and continuous darkness during winter days.
Example:
Compute the time of sunrise on 25 June at Wayne, New Jersey.
Latitude: \(40^{\circ}\!.9 \text{ North} \quad \phi=+40^{\circ}\!.9 \quad \sin \phi=+0.65474 \quad \cos \phi=+0.75585\)
Longitude: \(74^{\circ}\!.3 \text{ West} \quad \lambda=-74^{\circ}\!.3/15 = -4^{h}\!.953\)
For sunrise: \(z=90^{\circ} 50^{\prime} \quad \cos z = -0.01454\)
Sunrise occurs at \(9^{h} 26^{m}\) UT = \(5^{h} 26^{m}\) EDT
Alternative 3¶
The official citation:
United States Naval Observatory. Nautical Almanac Office. (19801991). Almanac for computers. Washington, D.C.: Nautical Almanac Office, United States Naval Observatory.
This has a section, Equation of Time and Time of Solar Transit. This computes local mean time (LMT) of noon.
We can offset this using the hour angle to compute sunrise and sunset. The hour angle requires the zenith, \(z\), declension of the sun, \(\delta\), and the latitude, \(\phi\). This amounts to the Alternative 2 solution, and doesn’t offer any benefit.
Equation of Time and Time of Solar Transit¶
This starts on Page B8.
The equation of time \(EqT\) is the hour angle of the true Sun minus the hour angle of the mean sun. Thus it is the difference: apparent solar (sundial) time minus mean solar (clock) time.
- \(N\)
Integer number of days since 0 January, \(0^{h}\) UT.
- \(t\)
Time since 0 January, \(0^{h}\) UT, in fractional days.
Approximation 1.
For higher accurancy, here is approximation 2.
In eq. (3), \(EqT\), the arctangent should yield a result in degrees that is in the same quadrant as \(\theta\). Near the end of the year \(\theta\) becomes greater than 360°. When this occurs the arctangent in eq. (3) should also be greater than 360°.
… [C]ompute \(EqT\) for \(t=N + \frac{12^{h}-\lambda}{24}\), where \(N\) is the day of the year … and \(\lambda\) is the longitude (east positive, west negative) expressed in hours. Then the local mean time (LMT) of transit is given to an accurace of \(\pm 2\) seconds by \(LMT = 12^{h} - EqT\). The univeral time of transit is then obtained with \(UT = LMT - \lambda\).
Procedure:
Compute \(t\) from \(N\) and \(\lambda\). \(t=N + \frac{12^{h}-\lambda}{24}\)
Compute \(EqT\) using one of the two approximations.
Compute LMT. \(LMT = 12^{h} - EqT\).
If needed, compute UT. \(UT = LMT - \lambda\).
Example:
Compute the time of solar transit at longitude \(73^{\circ}58^{\prime}\) West on 17 June 1990.