Sunrise equation

Last updated
A contour plot of the hours of daylight as a function of latitude and day of the year, using the most accurate models described in this article. It can be seen that the area of constant day and constant night reach up to the polar circles (here labeled "Anta. c." and "Arct. c."), which is a consequence of the earth's inclination. Hours of daylight vs latitude vs day of year with tropical and polar circles.svg
A contour plot of the hours of daylight as a function of latitude and day of the year, using the most accurate models described in this article. It can be seen that the area of constant day and constant night reach up to the polar circles (here labeled "Anta. c." and "Arct. c."), which is a consequence of the earth's inclination.
A plot of hours of daylight as a function of the date for changing latitudes. This plot was created using the simple sunrise equation, approximating the sun as a single point and does not take into account effects caused by the atmosphere or the diameter of the Sun.

The sunrise equation or sunset equation can be used to derive the time of sunrise or sunset for any solar declination and latitude in terms of local solar time when sunrise and sunset actually occur.

Contents

Formulation

It is formulated as:

where:

is the solar hour angle at either sunrise (when negative value is taken) or sunset (when positive value is taken);
is the latitude of the observer on the Earth;
is the sun declination.

Principles

The Earth rotates at an angular velocity of 15°/hour. Therefore, the expression , where is in degree, gives the interval of time in hours from sunrise to local solar noon or from local solar noon to sunset.

The sign convention is typically that the observer latitude is 0 at the equator, positive for the Northern Hemisphere and negative for the Southern Hemisphere, and the solar declination is 0 at the vernal and autumnal equinoxes when the sun is exactly above the equator, positive during the Northern Hemisphere summer and negative during the Northern Hemisphere winter.

The expression above is always applicable for latitudes between the Arctic Circle and Antarctic Circle. North of the Arctic Circle or south of the Antarctic Circle, there is at least one day of the year with no sunrise or sunset. Formally, there is a sunrise or sunset when during the Northern Hemisphere summer, and when during the Southern Hemisphere winter. For locations outside these latitudes, it is either 24-hour daytime or 24-hour nighttime.

Expressions for the solar hour angle

In the equation given at the beginning, the cosine function on the left side gives results in the range [-1, 1], but the value of the expression on the right side is in the range . An applicable expression for in the format of Fortran 90 is as follows:

omegao = acos(max(min(-tan(delta*rpd)*tan(phi*rpd), 1.0), -1.0))*dpr

where omegao is in degree, delta is in degree, phi is in degree, rpd is equal to , and dpr is equal to .

The above expression gives results in degree in the range . When , it means it is polar night, or 0-hour daylight; when , it means it is polar day, or 24-hour daylight.

Hemispheric relation

Suppose is a given latitude in Northern Hemisphere, and is the corresponding sunrise hour angle that has a negative value, and similarly, is the same latitude but in Southern Hemisphere, which means , and is the corresponding sunrise hour angle, then it is apparent that

,

which means

.

The above relation implies that on the same day, the lengths of daytime from sunrise to sunset at and sum to 24 hours if , and this also applies to regions where polar days and polar nights occur. This further suggests that the global average of length of daytime on any given day is 12 hours without considering the effect of atmospheric refraction.

Generalized equation

Sextant sight reduction procedure showing solar altitude corrections for refraction and elevation. Corrections for Sextant Altitude.en.jpg
Sextant sight reduction procedure showing solar altitude corrections for refraction and elevation.

The equation above neglects the influence of atmospheric refraction (which lifts the solar disc — i.e. makes the solar disc appear higher in the sky — by approximately 0.6° when it is on the horizon) and the non-zero angle subtended by the solar disc — i.e. the apparent diameter of the sun — (about 0.5°). The times of the rising and the setting of the upper solar limb as given in astronomical almanacs correct for this by using the more general equation

with the altitude angle (a) of the center of the solar disc set to about −0.83° (or −50 arcminutes).

The above general equation can be also used for any other solar altitude. The NOAA provides additional approximate expressions for refraction corrections at these other altitudes. [1] There are also alternative formulations, such as a non-piecewise expression by G.G. Bennett used in the U.S. Naval Observatory's "Vector Astronomy Software". [2]

Complete calculation on Earth

The generalized equation relies on a number of other variables which need to be calculated before it can itself be calculated. These equations have the solar-earth constants substituted with angular constants expressed in degrees.

Calculate current Julian day

where:

is the number of days since Jan 1st, 2000 12:00.
is the Julian date;
2451545.0 is the equivalent Julian year of Julian days for Jan-01-2000, 12:00:00.
0.0008 is the fractional Julian Day for leap seconds and terrestrial time (TT).
TT was set to 32.184 sec lagging TAI on 1 January 1958. By 1972, when the leap second was introduced, 10 sec were added. By 1 January 2017, 27 more seconds were added coming to a total of 69.184 sec. 0.0008=69.184 / 86400 without DUT1.
The operation rounds up to the next integer day number n.

Mean solar time

where:

is an approximation of mean solar time at integer expressed as a Julian day with the day fraction.
is the longitude (west is negative, east is positive) of the observer on the Earth;

Solar mean anomaly

where:

M is the solar mean anomaly used in the next three equations.

Equation of the center

where:

C is the Equation of the center value needed to calculate lambda (see next equation).
1.9148 is the coefficient of the Equation of the Center for the planet the observer is on (in this case, Earth)

Ecliptic longitude

where:

λ is the ecliptic longitude.
102.9372 is a value for the argument of perihelion.

Solar transit

where:

Jtransit is the Julian date for the local true solar transit (or solar noon).
2451545.0 is noon of the equivalent Julian year reference.
is a simplified version of the equation of time. The coefficients are fractional days.

Declination of the Sun

where:

is the declination of the sun.
23.4397° is Earth's maximum axial tilt toward the sun [3]

Hour angle

This is the equation from above with corrections for atmospherical refraction and solar disc diameter.

where:

ωo is the hour angle from the observer's meridian;
is the north latitude of the observer (north is positive, south is negative) on the Earth.

For observations on a sea horizon needing an elevation-of-observer correction, add , or to the −0.833° in the numerator's sine term. This corrects for both apparent dip and terrestrial refraction. For example, for an observer at 10,000 feet, add (−115°/60) or about −1.92° to −0.833°. [4]

Calculate sunrise and sunset

where:

Jrise is the actual Julian date of sunrise;
Jset is the actual Julian date of sunset.

Example of implementation in Python

#!/usr/bin/python3importloggingfromdatetimeimportdatetime,timedelta,timezone,tzinfofrommathimportacos,asin,ceil,cos,degrees,fmod,radians,sin,sqrtfromtimeimporttimelog=logging.getLogger()def_ts2human(ts:int|float,debugtz:tzinfo|None)->str:returnstr(datetime.fromtimestamp(ts,debugtz))defj2ts(j:float|int)->float:return(j-2440587.5)*86400defts2j(ts:float|int)->float:returnts/86400.0+2440587.5def_j2human(j:float|int,debugtz:tzinfo|None)->str:ts=j2ts(j)returnf'{ts} = {_ts2human(ts,debugtz)}'def_deg2human(deg:float|int)->str:x=int(deg*3600.0)num=f'∠{deg:.3f}°'rad=f'∠{radians(deg):.3f}rad'human=f'∠{x//3600}°{x//60%60}{x%60}″'returnf'{rad} = {human} = {num}'defcalc(current_timestamp:float,f:float,l_w:float,elevation:float=0.0,*,debugtz:tzinfo|None=None,)->tuple[float,float,None]|tuple[None,None,bool]:log.debug(f'Latitude               f       = {_deg2human(f)}')log.debug(f'Longitude              l_w     = {_deg2human(l_w)}')log.debug(f'Now                    ts      = {_ts2human(current_timestamp,debugtz)}')J_date=ts2j(current_timestamp)log.debug(f'Julian date            j_date  = {J_date:.3f} days')# Julian day# TODO: ceil ?n=ceil(J_date-(2451545.0+0.0009)+69.184/86400.0)log.debug(f'Julian day             n       = {n:.3f} days')# Mean solar timeJ_=n+0.0009-l_w/360.0log.debug(f'Mean solar time        J_      = {J_:.9f} days')# Solar mean anomaly# M_degrees = 357.5291 + 0.98560028 * J_  # Same, but looks uglyM_degrees=fmod(357.5291+0.98560028*J_,360)M_radians=radians(M_degrees)log.debug(f'Solar mean anomaly     M       = {_deg2human(M_degrees)}')# Equation of the centerC_degrees=1.9148*sin(M_radians)+0.02*sin(2*M_radians)+0.0003*sin(3*M_radians)# The difference for final program result is few milliseconds# https://www.astrouw.edu.pl/~jskowron/pracownia/praca/sunspot_answerbook_expl/expl-4.html# e = 0.01671# C_degrees = \#     degrees(2 * e - (1 / 4) * e ** 3 + (5 / 96) * e ** 5) * sin(M_radians) \#     + degrees(5 / 4 * e ** 2 - (11 / 24) * e ** 4 + (17 / 192) * e ** 6) * sin(2 * M_radians) \#     + degrees(13 / 12 * e ** 3 - (43 / 64) * e ** 5) * sin(3 * M_radians) \#     + degrees((103 / 96) * e ** 4 - (451 / 480) * e ** 6) * sin(4 * M_radians) \#     + degrees((1097 / 960) * e ** 5) * sin(5 * M_radians) \#     + degrees((1223 / 960) * e ** 6) * sin(6 * M_radians)log.debug(f'Equation of the center C       = {_deg2human(C_degrees)}')# Ecliptic longitude# L_degrees = M_degrees + C_degrees + 180.0 + 102.9372  # Same, but looks uglyL_degrees=fmod(M_degrees+C_degrees+180.0+102.9372,360)log.debug(f'Ecliptic longitude     L       = {_deg2human(L_degrees)}')Lambda_radians=radians(L_degrees)# Solar transit (julian date)J_transit=2451545.0+J_+0.0053*sin(M_radians)-0.0069*sin(2*Lambda_radians)log.debug(f'Solar transit time     J_trans = {_j2human(J_transit,debugtz)}')# Declination of the Sunsin_d=sin(Lambda_radians)*sin(radians(23.4397))# cos_d = sqrt(1-sin_d**2) # exactly the same precision, but 1.5 times slowercos_d=cos(asin(sin_d))# Hour anglesome_cos=(sin(radians(-0.833-2.076*sqrt(elevation)/60.0))-sin(radians(f))*sin_d)/(cos(radians(f))*cos_d)try:w0_radians=acos(some_cos)exceptValueError:returnNone,None,some_cos>0.0w0_degrees=degrees(w0_radians)# 0...180log.debug(f'Hour angle             w0      = {_deg2human(w0_degrees)}')j_rise=J_transit-w0_degrees/360j_set=J_transit+w0_degrees/360log.debug(f'Sunrise                j_rise  = {_j2human(j_rise,debugtz)}')log.debug(f'Sunset                 j_set   = {_j2human(j_set,debugtz)}')log.debug(f'Day length                       {w0_degrees/(180/24):.3f} hours')returnj2ts(j_rise),j2ts(j_set),Nonedefmain():logging.basicConfig(level=logging.DEBUG)latitude=33.00801longitude=35.08794elevation=0print(calc(time(),latitude,longitude,elevation,debugtz=timezone(timedelta(hours=3),'fake-zone')))if__name__=='__main__':main()

See also

Related Research Articles

<span class="mw-page-title-main">Gyrocompass</span> Type of non-magnetic compass based on the rotation of the Earth

A gyrocompass is a type of non-magnetic compass which is based on a fast-spinning disc and the rotation of the Earth to find geographical direction automatically. A gyrocompass makes use of one of the seven fundamental ways to determine the heading of a vehicle. A gyroscope is an essential component of a gyrocompass, but they are different devices; a gyrocompass is built to use the effect of gyroscopic precession, which is a distinctive aspect of the general gyroscopic effect. Gyrocompasses, such as the fibre optic gyrocompass are widely used to provide a heading for navigation on ships. This is because they have two significant advantages over magnetic compasses:

<span class="mw-page-title-main">Astronomical coordinate systems</span> System for specifying positions of celestial objects

In astronomy, coordinate systems are used for specifying positions of celestial objects relative to a given reference frame, based on physical reference points available to a situated observer. Coordinate systems in astronomy can specify an object's position in three-dimensional space or plot merely its direction on a celestial sphere, if the object's distance is unknown or trivial.

In geometry, a solid angle is a measure of the amount of the field of view from some particular point that a given object covers. That is, it is a measure of how large the object appears to an observer looking from that point. The point from which the object is viewed is called the apex of the solid angle, and the object is said to subtend its solid angle at that point.

In control theory and signal processing, a linear, time-invariant system is said to be minimum-phase if the system and its inverse are causal and stable.

In physics, the Hamilton–Jacobi equation, named after William Rowan Hamilton and Carl Gustav Jacob Jacobi, is an alternative formulation of classical mechanics, equivalent to other formulations such as Newton's laws of motion, Lagrangian mechanics and Hamiltonian mechanics.

The solar zenith angle is the zenith angle of the sun, i.e., the angle between the sun’s rays and the vertical direction. It is the complement to the solar altitude or solar elevation, which is the altitude angle or elevation angle between the sun’s rays and a horizontal plane. At solar noon, the zenith angle is at a minimum and is equal to latitude minus solar declination angle. This is the basis by which ancient mariners navigated the oceans.

The solar azimuth angle is the azimuth of the Sun's position. This horizontal coordinate defines the Sun's relative direction along the local horizon, whereas the solar zenith angle defines the Sun's apparent altitude.

<span class="mw-page-title-main">Phasor</span> Complex number representing a particular sine wave

In physics and engineering, a phasor is a complex number representing a sinusoidal function whose amplitude, and initial phase are time-invariant and whose angular frequency is fixed. It is related to a more general concept called analytic representation, which decomposes a sinusoid into the product of a complex constant and a factor depending on time and frequency. The complex constant, which depends on amplitude and phase, is known as a phasor, or complex amplitude, and sinor or even complexor.

A resistor–inductor circuit, or RL filter or RL network, is an electric circuit composed of resistors and inductors driven by a voltage or current source. A first-order RL circuit is composed of one resistor and one inductor, either in series driven by a voltage source or in parallel driven by a current source. It is one of the simplest analogue infinite impulse response electronic filters.

In mathematics, a change of variables is a basic technique used to simplify problems in which the original variables are replaced with functions of other variables. The intent is that when expressed in new variables, the problem may become simpler, or equivalent to a better understood problem.

In the mathematical description of general relativity, the Boyer–Lindquist coordinates are a generalization of the coordinates used for the metric of a Schwarzschild black hole that can be used to express the metric of a Kerr black hole.

The Rabi problem concerns the response of an atom to an applied harmonic electric field, with an applied frequency very close to the atom's natural frequency. It provides a simple and generally solvable example of light–atom interactions and is named after Isidor Isaac Rabi.

<span class="mw-page-title-main">Voigt effect</span>

The Voigt effect is a magneto-optical phenomenon which rotates and elliptizes linearly polarised light sent into an optically active medium. The effect is named after the German scientist Woldemar Voigt who discovered it in vapors. Unlike many other magneto-optical effects such as the Kerr or Faraday effect which are linearly proportional to the magnetization, the Voigt effect is proportional to the square of the magnetization and can be seen experimentally at normal incidence. There are also other denominations for this effect, used interchangeably in the modern scientific literature: the Cotton–Mouton effect and magnetic-linear birefringence, with the latter reflecting the physical meaning of the effect.

<span class="mw-page-title-main">Salah times</span> Timing of Islamic prayers

Salat times are prayer times when Muslims perform salat. The term is primarily used for the five daily prayers including the Friday prayer, which takes the place of the Dhuhr prayer and must be performed in a group of aibadat. Muslims believe the salah times were revealed by Allah to Muhammad.

In the differential geometry of surfaces, a Darboux frame is a natural moving frame constructed on a surface. It is the analog of the Frenet–Serret frame as applied to surface geometry. A Darboux frame exists at any non-umbilic point of a surface embedded in Euclidean space. It is named after French mathematician Jean Gaston Darboux.

In mathematics, vector spherical harmonics (VSH) are an extension of the scalar spherical harmonics for use with vector fields. The components of the VSH are complex-valued functions expressed in the spherical coordinate basis vectors.

<span class="mw-page-title-main">Geographical distance</span> Distance measured along the surface of the Earth

Geographical distance or geodetic distance is the distance measured along the surface of the Earth, or the shortest arch length.

<span class="mw-page-title-main">Capstan equation</span> Relates the hold-force to the load-force if a flexible line is wound around a cylinder

The capstan equation or belt friction equation, also known as Euler-Eytelwein formula, relates the hold-force to the load-force if a flexible line is wound around a cylinder.

An electric dipole transition is the dominant effect of an interaction of an electron in an atom with the electromagnetic field.

The Bowring series of the transverse mercator published in 1989 by Bernard Russel Bowring gave formulas for the Transverse Mercator that are simpler to program but retain millimeter accuracy.

References

  1. NOAA (U.S. Department of Commerce). "Solar Calculation Details". ESRL Global Monitoring Laboratory - Global Radiation and Aerosols.
  2. "Correction Tables for Sextant Altitude". www.siranah.de.
  3. "Earth Fact Sheet".
  4. The exact source of these numbers are hard to track down, but Notes on the Dip of the Horizon provides a description yielding one less significant figure, with another page in the series providing -2.075.