[Physics] How to calculate how many degrees the Sun is from the horizon

suntime

Given a position on earth (latitude/longitude, and maybe also height relative to sea-level) and a time, what's the algorithm to find how many degrees the Sun is from the horizon?

And the same for the opposite, given a location and date, how do you calculate the time when the Sun will be such-and-such degrees from the horizon (such as to calculate astronomical dusk and dawn)?

Best Answer

There are a number of steps involved, each of which are straightforward individually, but when combined make the process a bit tedious to do by hand. The following steps are only a conceptual outline, and ignore a number of important considerations, including but not limited to the ellipticity of the Earth's orbit. The outline of the algorithm is:

Calculate the position of the Sun in ecliptic coordinates. By definition, the ecliptic latitude of the Sun is 0. If you are willing to accept the inaccuracy introduced by pretending the Earths orbit is circular rather than elliptical, and can look up the ecliptic coordinates of the Sun for some time relatively close to the time for which you want the position, you can make the approximation $$\lambda = \lambda_D + \frac{360^\circ}{365.242191\, \mathrm{days}} \Delta D$$ where $\lambda_D$ is the ecliptic longtitude of the sun at date $D$, and $\Delta D$ is the number of days between the time you want and date $D$.
Convert the ecliptic coordinates to equatorial coordinates. Because the ecliptic latitude of the Sun is 0, this reduces to: $$\alpha = \tan^{-1} \left[\frac{\sin \lambda \cos \epsilon}{\cos \lambda}\right]$$ $$\delta = \sin^{-1} [\sin \epsilon \sin \lambda]$$ where $\epsilon$ is the obliquity of the ecliptic, currently about 23.4378 degrees, but it does change over time.
Now you need to convert the right ascension in ecliptic coordinates to an hour angle at the time and longitude of interest. To do this, you need to calculate the local sidereal time. (You can do this by first calculating the Greenwich Sidereal Time (GST), and converting to local sidereal time. Sidereal time and solar time agree at the autumn equinox, so look up the autumn equinox in UT, and calculate the time between then and the time of interest. A sidereal year is one day longer than a solar year, so convert that duration from solar time to sidereal time by multiplying by $\frac{365.242191+1}{365.242191}$. Then convert Greenwich sidereal time to local sidereal time by subtracting the (west) longitude. Be careful with the units here: 24 hours is 360 degrees is $2 \pi$ radians.) Then you can calculate the hour angle: $$H=\mathrm{LST} - \alpha$$.
Once you have the hour angle, $H$, you can calculate the altitude, $a$ (which is half of the conversion to horizon coordinates): $$a=\sin^{-1} [\sin \delta \sin \phi + \cos \delta \cos \phi \cos H]$$ where is $\phi$ is the geographic latitude of interest.

For a treatment sufficient for actual calculation, you can start with Practical Astronomy with your Calculator or Spreadsheet by Duffett-Smith & Zwart. For higher precision, you can go to Astronomical Algorithms by Meeus.

Related Solutions

[Physics] actually meant by ‘sun set’ and ‘sun rise’ times, when taking into account the mirage due to light bending in the atmosphere

See http://en.wikipedia.org/wiki/Sunset. The sunset is defined when the trailing edge of the sun disappears below the horizon, even though refraction means the sun actually went below the horizon some minutes before. Likewise sunrise is defined as the point when the leading edge of the sun is first visible even though this is some minutes before the sun actually rises above the horizon.

You ask where the sun really is: when the sun sets and rises it is about one solar diameter below where it appears to be, i.e. about half a degree. See http://en.wikipedia.org/wiki/Atmospheric_refraction for the details.

[Physics] Reverse Sun position algorithm

https://stackoverflow.com/questions/257717/position-of-the-sun-given-time-of-day-and-lat-long has a good set of descriptions of how to do approximately what you want. In short, you either want to use a code or look up the sun's ephemeris.

Using pyephem, for example:

import ephem
o = ephem.Observer()
# lon/lat for LA
o.lon = "118:15:00"    
o.lat = "34:03"
sun = ephem.Sun(o)
sun.alt
sun.az

But since you want the inverse, you can get the alt/az as a function of time:

import datetime
dt5 = datetime.timedelta(minutes=5)
time = datetime.datetime.now()
alt,az,times = [],[],[]
for ii in xrange(100): # do 100 steps of 5 minutes
    time = time+dt5
    o.date = time.strftime("%Y/%m/%d %H:%M:%S")
    sun.compute(o)
    alt.append(sun.alt)
    az.append(sun.az)
    times.append(o.date)

You can then invert this data to find time as a function of alt/az. That doesn't quite answer your question, though - I'll keep an eye out for a direct answer.

Edit: A direct answer.

The Hour Angle HA is the current time referenced to Solar noon. It can be computed by:

$$\cos(HA) = \frac{\sin(alt_\odot) - \sin(\delta_\odot) \sin(Latitude)}{cos(\delta_\odot)} $$ Where $alt_\odot$ is the altitude or elevation angle of the sun and $\delta_\odot$ is the declination of the sun and can be approximated by: $$\delta_\odot = -23.44^o \cos\left(\frac{360^o}{365} (N+10)\right)$$ where $N$ is the ordinal date and is 1 at January 1st

Best Answer

Related Solutions

[Physics] actually meant by ‘sun set’ and ‘sun rise’ times, when taking into account the mirage due to light bending in the atmosphere

[Physics] Reverse Sun position algorithm

Related Question