|
Algorithms and Formulas |
How to find the times of sunrise
and sunset,
by Vanni
Tomezzoli, March 17, 2001.
The sunrise occurs when the upper limb of the Sun disc is visible at the horizon, towards east, at a location whose elevation is reduced to the sea level - while the sunset occurs in the same circumstances, but in the opposite direction, towards the western horizon.
Top of pageThere is not a direct way to compute the times of sunrise and sunset for any location, at a given date; this can be performed either by interpolation or - as in the case of Interactive Sunriset - by iteraction. This is due to one of the variables involved: time, that is the solution itself. The calculation is faster with the interpolation, expecially when requested for several days, while the iteraction should be more accurate. Anyway, an absolute accuracy cannot exist, as the apparent position of the Sun is affected by the local atmospheric refraction, whose effects cannot be predicted.
The information required for the computation of the times of sunrise and sunset are:
Next, you must know the orbital elements to find out the position of the Sun at the requested date. Like any celestial chart, which reports the date (“epoch”) of the mapped sky, the astronomical calculation refers to a particular moment. For example, “Interactive Sunriset” refers to the position of the Sun at epoch 2000.0. This is also called the standard equinox 2000.0, which refers to the position of the vernal equinox at 12 hours (noon) of the Greenwich meridian, on January 1, 2000. (The zero digit specified after the decimal point indicates the beginning of the year, therefore “epoch 2000.0” stands for the beginning of the year 2000). The orbital elements are relative to this date.
The equinox is a straight line where the plane of the ecliptic (the orbit of the Earth around the Sun) meets the equatorial plane of the Earth. On March 21 the Sun is placed on this line, where the ecliptic arise towards the northern hemisphere (ascending node), while on September 23 the Sun is placed on the same line, but in the opposite direction. The equinox itself is moving slightly with a cycle of about 26,000 years and is perpendicular to the Earth's axis of rotation. This motion is named the “precession of the equinoxes”.
The Universal Time (UTC) is the mean time at the Greenwich meridian. The calculation needs the value of different elements at the time of the referred equinox at epoch 2000.0, and the quantities involved to upgrade the same values to the requested time (e.g., today). The calculation referred to any epoch are valid for several centuries after and before the same epoch, it is determinant therefore to know the time lapsed from the equinox of the epoch of reference to the date you request.
You should refer to the time using the Julian Day (JD), a time system used
for astronomical purposes.
The Julian Day must not be confused with the Julian calendar. Time in JD
is not divided in years, months and days of the month
as in our civil life - but only Julian days
and a fractional part, represents hours of the day, and they are referred to the
Greenwich Mean Time, with the day beginning at noon instead of midnight.
When you connect to the home page of my web site you can see the Julian day
relative to the moment of your connection, following the date and UTC.
Subtracting the JD of the epoch of reference from the JD of the requested date,
(e.g., the JD for today), you obtain the number of days lapsed since the epoch
of the equinox.
Top of page
T Dynamic Time
The time as shown by our clocks is determined by the astronomers, who report after their observations, measuring the Earth rotation on its axis. The time of rotation of the Earth is not regular, therefore it isn't a constant value. Only atomic clocks are able to report such irregularities, which can be quantified on about one second per year. This shift isn't predictable, and the value is available after the observation. In astronomical calculation you need to refer to a constant interval of time, called Dynamic Time, which must be adjusted to the Universal Time (UTC), which refers to the effective rotation of the Earth and affects the synchronization of our clocks. For the year 2001 the predicted value for the Dynamyc Time is about 65 seconds, which means that the constant Dynamic Time if 65 seconds in advance of the actual UTC.
T = DT - UTC = 65 seconds
for the year 2001.
Find the Julian Day for the requested date (e.g. today) at 0 hours UTC (that is, the midnight at Greenwich). Store this value as you'll need it later, at the end of the iteration, and refer to this as the starting JD (or JD of reference). Get the Julian Day from any astronomical almanac. If you want to determine the JD writing the software yourself, here is the algorithm effective for the Gregorian Calendar. Supply the day (a fractional part for the hours may be included, but this isn't requested at this moment), the month (1 to 12) and the year (four digits).
If the month is lesser than 3, add 12 to the number of the month and subtract 1 from the year.
n = int (year / 100)
JD = int (365.25 year) + int (30.6001 (month + 1)) + day + 1720996.5 - n + int (n / 4)
“int” refers to the integer part of the expressions enclosed within parenthesis, rounding to the lower integer value.
One Julian century includes 36,525 Julian days. Compute Julian centuries lapsed from the epoch 2000.0 to the date of reference (for example, a value of 0.25 for the variable T refers to the beginning of the year 2025; T is negative before the year 2000).
|
JD0 - 2451545
|
||
| T = |
|
|
|
36525
|
Compute the square and cube of T necessary for the secular terms: T2 and T3.
The sidereal time is synchronized to the rotation of the Earth around its axis, and the duration of one sidereal day is about 23 hours and 56 minutes. During this time the Earth proceeds for almost one degree in its revolution around the Sun, and must spend four minutes further to see the Sun at the same position of the previous day. By the other hand, the mean time, used for our civil life, is synchronized to the transit of the Sun at the meridian line. One solar day may be defined as the time lapsed between two consecutive middays (about 24 hours).
GST0 = 6.6973745583 + 2400.0513369072 T + 0.000025862 T2 - 0.00000000172 T3
Having the value for JD of reference found earlier, the expected value of JD for the sunrise may be adjusted as follows:
|
L
|
||
|
JDe = JD0 + 0.25 +
|
|
|
|
360°
|
while for the sunset may be adjusted as follows:
|
L
|
||
|
JDe = JD0 + 0.75 +
|
|
|
|
360°
|
These dummy values are not effective for the true time of sunrise and sunset,
but we first need them as the expected value for the necessary interation.
Referring to one of these values (the expected JD for the sunrise
or the sunset), the iteraction begins here.
Top of page
First you must know the Julian centuries lapsed from the standard equinox to the expected (dummy) JD, including a fractional part for the hours and minutes.
Compute a new value for T (Julian centuries) adapted for the expected time lapsed from the epoch 2000.0:
|
JDe + DTime - 2451545
|
||
| T = |
|
|
|
36525
|
You need also the new values for T2 and T3, in the same way we did for the JD of reference.
Next, find out the position of the Sun as follows:
L0 = 280.46646 + 36000.76983 T + 0.0003032 T2
The anomaly is the angular difference between a mean circular orbit and the true elliptic orbit.
M = 357.52911 + 35999.05029 T - 0.0001537 T2
C = (1.914602 - 0.004817 T - .000014 T2) sin M + (0.019993 - 0.000101 T) sin (2 M) + 0.000289 sin (3 M)
True longitude of the Sun
= L0 + C
v = M + C
The eccentricity is the ratio between the semi-major axis and the difference between the semi-major and semi-minor axis of the elliptic orbit of the Earth around the Sun.
e = 0.016708634 - 0.000042037 T - 0.0000001267 T2
The distance from a planet (here is the Earth) to the Sun, expressed in AU. AU (Astronomical Unit) is the mean distance of the Earth from the Sun, about 150 millions Kms.
|
1.000001018 (1 - e 2)
|
||
| R = |
|
|
|
1 + e cos v
|
|
0.266563889
|
||
| SD = |
|
|
|
R
|
Longitude of the Moon's ascending node = 125.04452 - 1934.136261 T + 0.0020708 T2 + T3 / 450000
Apparent longitude of the Sun (ecliptic) The ecliptic longitude of the Sun is corrected for the nutation and the aberration. The nutation is the deviation of the Earth's axis of rotation, referred to the precession of the equinox.
= - .00569 - .00478 sin
Obliquity of the eclipticrelative to the Earth's equator.
= 23.4392911 - 0.01300416667 T - 0.00000016389 T2 + 0.00000050361 T3
+ 0.00255625 cos
The ecliplic coordinate system used so far refers to the ecliptic plane, and must converted to the equatorial coordinates system; this is a projection of the Earth coordinates to the celestial sphere, and is synchronized with the rotation of the Earth.
The declination may be considered as the latitude (the angular distance of a celestial body from the celestial equator), while the right ascension is the longitude (the angular distance from the vernal equinox, expressed either in hours or degrees).
Declination of the Sunsin = sin sin
Right ascension of the Sun (in degrees)
y = cos sin
x = cos
tan = y / x
The value obtained here lies in the range from -90 to +90 degrees, but you need a value included in the range from 0 to 360 degrees. Adjust the result for the correct quadrant as follows:
if x is negative add 180° to the right ascension (second and third quadrant) otherwise if y is negative add 360° (fourth quadrant).
The angular distance in longitude from a celestial body to the local meridian line. May be expressed in degrees (0° to 360°) or hours (0 to 24 hours). Here we are considering in degrees.
The meridian is the line crossing the celestial sphere above an observer's location, from the North celestial pole and passing through the zenith (the highest point in the sky) and falling to the South point of the horizon, in the northern hemisphere; from the South celestial pole, passing through the zenith and falling to the North point of the horizon, in the southern hemisphere.
y = -sin (0°34' + SD) - sin 
sin
x = cos cos
O°34' is the mean value for the atmospheric refraction, a phenomenon which shifts the apparent position of celestial bodies, towards the zenith point. This effect is increased at the horizon. The quantity of refraction may change depending on the actual air pressure and temperature. An average amount of 0°34' is usually adopted (some people consider a value of 0°36'). Due to the effects of the refraction, the time of sunrise and sunset cannot be foreseen with absolute accuracy.
If you wish, you can add to 0°34' the result of the following formula:
0.03
where h is the elevation in meters from the level of the sea, but usually the times of sunrise and sunset are reduced for sea level.
In the following circumstances the Sun doesn't rise or set:
-1 the Sun is above the horizon (midnight's Sun).
+1 the Sun is below the horizon.
Otherwise,
For sunset:
cos H = y / x
For sunrise:
cos H = -(y / x)
at the time of the local sunrise or sunset.
|
H + L +
|
||
| GST = |
|
|
|
15
|
at the time of the local sunrise or sunset.
UTC = 0.997269566329875 (GST - GST0)
Adjust the Julian Day to the local mean time:
|
UTC
|
L
|
|||
| JD = JD0 + |
|
+
|
|
|
|
24
|
360°
|
This value must be compared to the expected (dummy) value:
JD - JDe
If the absolute result is higher than the requested accuracy, assume that the UTC found so far as the dummy time (JDe) for a new iteration; otherwise the UTC is the time of sunrise or sunset, you have got the solution and the loop ends here.
A value of 0.0001 Julian Days returns a level of accuracy better
than nine seconds - but sometimes lesser than the effects of the
local atmospheric refraction.
Top of page
Just convert UTC to local time, adding DST where applicable. The DST (Daylight Saving Time) is the time applied in most countries during the summer, when the clocks are shifted one hour forward.
The angular position of an object or a celestial body, projected to the horizon line, measured clockwise from the north, as follows:
Here the azimuth is the direction of the Sun at the horizon, at the moment of the sunrise or sunset.
|
sin |H|
|
||
| tan A = |
|
|
|
cos |H| sin L - tan
cos L
|
|H| is the absolute value of H (hour angle).
If azimuth is positive, subtract azimuth from 180°, otherwise consider the absolute value of the azimuth (that is, change sign).
At last, for the sunset, subtract azimuth from 360°.
Top of pageSome useful routines from “Interactive Sunriset”, follow. I wrote them in JavaScript; they must be properly converted to other computer programming languages, if necessary.
JavaScript and other computer languages, such as Microsoft Visual Basic for Windows, perform trigonometric functions in radians, therefore all angular values must be converted from degrees to radians or vice versa.
// conversion from degrees to radians
Rad = new Function
("Value", "return Value * Math.PI / 180");
// conversion from radians to degrees
Deg = new Function
("Value", "return Value * 180 / Math.PI");
Math.PI is a JavaScript internal function which returns the value of Greek Pi (3.141592...).
Example:
Declination = Deg(Math.asin(Math.sin(Rad(Obliquity)) * Math.sin(Rad(ALongitude))));
// reduction to a range (24 hours or 360 degrees)
Range = new Function
("Limit", "Degrees", "return
(Degrees / Limit - Math.floor(Degrees / Limit))
* Limit");
Example:
MLongitude = Range (360, 280.46646 + 36000.76983 * Centuries + .0003032 * Centr2);
The degrees in latitude an longitude are usually expressed in degrees and minutes of arc, but in calculations you must use decimal values for the fractional part.
/* Conversion to decimal from degrees and minutes of arc */
function Dec (Degrees, Minutes) { var sign = (Degrees < 0) ? -1 : 1; return Math.floor (Degrees) + sign * (Minutes / 60); }
Here is the opposite operation, useful to convert the fractional part of the results
to hours and
minutes, for example, the time of sunrise and sunset.
/*
Conversion from decimal to degrees
and minutes of arc
(or hours and minutes)
*/
var Degs = 0; // degrees (integer part)
var Mins = 0; // minutes of arc
function Dms (Value) {
Degs = Math.floor(Value);
var temp = (Math.abs (Value - Degs)) * 60;
Mins = Math.floor(temp + .5);
if (Mins == 60) {
Mins = 0;
++Degs;
}
}
The year, month and day of the Gregorian Calendar must be supplied. The parameter day can be supplied with a fractional part for the hour of the day starting at midnight.
Example for December 25, 2000 at 18 hours UTC:
JD = JDay (2000, 12, 25.75);
Example for May 3, 2001 at 0 hours UTC:
JD = JDay (2001, 5, 3);
Here is the code for the function:
// Conversion to Julian Day
function JDay (year, month, day) {
if (month < 3) {
year -= 1;
month += 12;
}
var a = Math.floor (year / 100);
var b = 2 - a + Math.floor (a / 4);
return Math.floor (365.25 * year) +
Math.floor (30.6001 * (month + 1)) +
day + 1720994.5 + b;
}
Top of pageAstronomical Algorithms by Jean Meeus. Willmann-Bell, Inc., Second Edition, 1998. Most formulas were deducted from this essential book.
Textbook on Spherical Astronomy by W. M. Start, revised by R. M. Green. Cambridge University Press. Detailed information on the celestial mechanics.
The URL of Sunriset is:
© Copyright 2001 - Vanni Tomezzoli, Italy. All rights reserved.