Notes

airmass.org calculates the position of celestial objects as seen from any selected location on the surface of the Earth, over the course of both a night and a year. It also calculates the rise and set times of the Sun and the Moon. These notes describe how it works. Astronomical Algorithms by Jean Meeus is the basis for a lot of the calculations, and other sources of algorithms are noted where used.

airmass.org is developed in the UK, using Google Chrome running on Linux, and tested on Mozilla Firefox and mobile browsers. I'm not aware of any bugs right now but they could be present and affect other systems and (because of Javascript's localisation of dates) in other timezones. Please let me know if you encounter problems.

Object altitude

To calculate the altitude and azimuth of a celestial object from its Right Ascension and Declination, the local sidereal time is required, to put into the equations

sin(altitude)=sin(latitude)*sin(δ) + cos(latitude)*cos(δ)*cos(LHA)
cos(azimuth)=(sin(δ)-sin(altitude)*sin(latitude))/(cos(altitude)*cos(latitude))

where LHA is the local hour angle in degrees, obtained by subtracting the Right Ascension from the local sidereal time, and δ is the declination of the object. The local sidereal time is calculated using the method described at http://aa.usno.navy.mil/faq/docs/GAST.php (currently down; archived copy) neglecting the correction for nutation which amounts to only 1.1s at most.

The calculated altitude is corrected for atmospheric refraction, which means that objects reach an apparent zenith distance of 90° when their true zenith distance is already ≈90.5°. The equation used is from Sæmundsson 1986 as quoted by Meeus.

When an object's coordinates are resolved from its name, airmass.org also retrieves the proper motions in Right Ascension and Declination if known. It uses these to determine the RA and dec at the epoch of the observation, taken as the start of the year.

Precession is then calculated to convert the RA and dec (at epoch 2000.0 if the proper motion is negligible or not measured) from equinox 2000.0 to the equinox at Jan 1 of the year of the selected date. Proper motion typically has a negligible effect, even over many centuries. Precession causes a 50 arcsec shift per year in the coordinate system, which for J2000 coordinates observed in 2018 is 15 arcmin or 1 minute in RA.

To calculate the airmass, the simplest equation assumes a plane parallel atmosphere and is X=sec(z), where z is the zenith distance. airmass.org corrects for the curvature of the Earth using the equation of Rozenberg (1966), which is

X = (cos(z) + 0.25e-11cos(z))-1

Sunset and twilight times

To calculate sunset and sunrise times, airmass.org accounts for refraction and the elevation of the observing site by determining when the altitude of the centre of the Sun's disc is (-0.83 - (elevation)0.5/60) °, where elevation is in metres above sea level. The radius of the solar disc is 0.25° and refraction raises the apparent elevation by 0.58°. The altitude correction term implies that, at 3600m above sea level, the zenith distance of the horizon is increased by 1°.

Once the Sun is below the horizon, refraction no longer matters, and airmass.org does not account for elevation either, so the boundaries of civil, nautical and astronomical twilights are just calculated as the times when the altitude of the Sun reaches -6, -12 and -18°.

The equations for calculating the time at which the Sun reaches a given altitude are then:

MST = n - (longitude/360)
M = (357.5291 + 0.98560028*MST)%360
C=1.9148*sin(M) + 0.2*sin(2M) + 0.0003*sin(3M)
λ=(M+C+180+102.9372)%360
J_transit = jd2000 + 0.5 + MST + (0.0053*sin(M)) - (0.0069*sin(2λ))
sin(δ) = sin(λ) * sin(23.44)
cos(ω) = (sin(el) - sin(δ)*sin(latitude)) / (cos(latitude) * cos(δ))
J_event = J_transit ± ω/360

where

Everything on airmass.org is night-centric, so it calculates sunset and the starts of twilights using the relevant value of n for the date of the start of the night, while the ends of the twilights and the time of sunrise are calculated using n+1.

Moon phase and position

A really precise calculation of the Moon's position requires the calculation of hundreds of periodic terms to account for all the myriad gravitational forces acting on the Moon. Meeus gives a simplified set of calculations in chapter 45 of Astronomical Algorithms, with 59 periodic terms in longitude and 30 in latitude. For the sake of faster calculations, airmass.org simplifies this further by taking the first 13 periodic terms in the Moon's longitude and the first 10 in latitude.

The calculation of moonrise and moonset is done iteratively, as its RA and dec are changing constantly. airmass.org calculates the Moon's position every two minutes throughout the 24 hour period it is looking at, and from this, it finds the times at which its altitude changes from positive to negative or vice versa. It then linearly interpolates to estimate the time at which the altitude was zero.

Lunar calculations are made firstly in the geocentric frame and then corrected to the topocentric position. Refraction is then accounted for, and finally the rise and set times determined using the horizon corrected for the elevation of the observing location. Moon rise and set times should be accurate to within a few minutes.

The times of full moons throughout the year are calculated using the method described in Meeus's chapter 47, but using only the first 14 of 25 periodic correction terms, and omitting the planetary corrections. These simplifications introduce a maximum error of 2-3 minutes.

Many of the calculations rely on the Julian Date. All dates are stored internally as the number of milliseconds since 1 Jan 1970, 00:00 UT, so the Julian Date can be calculated by simply dividing that by 86400000 (the number of milliseconds in a day), and then adding 2440587.5, the Julian Date at Unix time zero. The difference between Universal Time and Coordinated Universal Time, which is 37 seconds as of July 2020, and increases at irregular intervals when leap seconds are added, is currently ignored in all calculations.

Parallactic angle

The parallactic angle is the angle between the meridian and the great circle passing through an object and the celestial poles. Aligning a spectrograph slit at the parallactic angle avoids differential refraction along the slit. airmass.org calculates the angle using the following equation, and displays it in the dynamic display on the daily chart as well as in the extended data for each object.

HA = RA - 15*LAST
PA = (360-arctan(sin(HA)/(tan(latitude)*cos(dec)-sin(dec)*cos(HA))))%360

where RA and dec are the celestial coordinates of the object, LAST is the local apparent sidereal time in hours, and HA is the hour angle.

Heliocentric velocity correction

The heliocentric velocity correction is given in the extended data. The algorithm is adapted from the rv command within starlink. The data shows you the total correction to apply, as well as the components due to the Earth's rotation and orbital velocity. The value given should be subtracted from a measured radial velocity to obtain a heliocentric radial velocity. The calculation should provide an accuracy of ±0.01km/s.

Solar system objects

Under "advanced options..." you can look up solar system objects from the IMCCE Virtual Observatory database. Positions are downloaded from the database for the current year only, and the RA and Dec are calculated at UTC midnight. The daily chart positions are interpolated using the motions on the sky returned with the ephemeris. For fast-moving objects close to Earth this introduces some imprecision but remains accurate enough to find out when the object is best seen. You can see this if you look up the Moon - the interpolated line will appear close to, but not exactly on, the calculated Moon line.

Benchmarking

I've compared the numbers from airmass.org to a few other almanac services. I randomly chose July 9 2018 and Paranal observatory, looking at NGC 5189 (link to calculation). Here is a table comparing results from airmass.org (calculated on 16 April 2021, internal version ID b93bc92), with results from the ING's staralt, Heavens Above, and ESO's skycalc.

observability.dateheavens aboveeso skycalcstaralt
Sunset18:15:3418:0818:1618:15
End of civil twilight18:31:4518:32
End of nautical twilight18:59:5819:0119:01
End of astronomical twilight19:27:4519:2819:2819:27
Start of astronomical twilight06:05:2006:0606:0606:05
Start of nautical twilight06:33:0606:3306:33
Start of civil twilight07:01:1707:01
Sunrise07:17:2707:2607:1807:18
Length of night (sunset-sunrise)13:0213:1813.013:03
Length of astronomical night10:3810:3810.610:38
LST at civil midnight18:3018:30
Moon altitude at midnight-58.9°
Distance from object124.5°123
Illuminated fraction0.130.140.1220.12
Days past full12
Sets15:17:2615:0115:09
Rises04:21:3004:3604:2804:40*
object max alt48.648*
airmass1.331.35*
parallactic angle at midnight8679
hours above 30 degrees in dark:4h22m
* = read from chart

There is generally fair agreement between the various pages. Heavens Above's discrepant sunset and sunrise times are because it does not correct for altitude, where the other options all do. All other solar phenomena agree to within a minute; actual variation in sunrise and sunset due to varying atmospheric conditions along the line of sight to the horizon is likely to be larger than this.

Lunar phenomena show much more variation. Moon rise and set times have a maximum spread of 14-18 minutes. Possible reasons for this include:

The heliocentric velocity correction calculated by airmass.org should be accurate to about ±100 m/s. For NGC 5189, at 12:00UT on 2018 July 09, it is +14.985 km/s. The barycentric correction for the same time as given by the algorithm of Wright & Eastman (2014), accurate to a few mm/s, is -14.990 km/s (Wright & Eastman use the opposite sign convention to rv; airmass.org adopts the latter).

Some further tests to check the accuracy of airmass.org's calculations can be run by going to the unit tests page. All calculations are done by your browser so the results may differ between systems, but should be very close to the benchmark values.

Accuracy

Most of the equations used in airmass.org cannot be used indefinitely, as they contain polynomial expressions in time which will diverge far from 2000. In case you should wish to use the site for dates very far into the future or past, the validity of the results is roughly as follows:

Celestial objects: the correction for precession is valid for at least 100 centuries. airmass.org covers dates from 0-9999AD (at least with the date support in Chrome) and so the precessed coordinates throughout that time should be accurate. The proper motion calculation is a simple method which will become inaccurate over long time spans for stars very close to the celestial poles.

Local sidereal time: airmass.org uses a method from the US Naval Observatory which they say loses precision at 0.1s/century. Thus, by 9999 (the last year for which the HTML5 date input works), the discrepancy should only be 8s.

Sun: Solar event calculations rely on calculations of the ecliptic longitude and latitude of the Sun. The equation for the latitude is accurate to within 1" for 1000BC-3000AD; longitude calculations are accurate to about 30" over a similar period.

Moon: Lunar calculations omit smaller periodic terms to keep the calculations fast. This limits the accuracy of lunar timings to a few minutes. Any secular trend should come from the secular variation of the Moon itself, an effect probably negligible within the 0-9999 date range available.

Light pollution

The darkness of the night sky at the selected observing location is taken from the World Atlas of the Night Sky Brightness, as updated in 2016. The amount of light pollution is shown in terms of the natural sky brightness: "1.0 × natural sky brightness" means that the night sky appears twice as bright as it does from the darkest sites on Earth. Site factors might make the night sky darker than the atlas suggests; for example, a cloud layer below the altitude of the observatory would make the sky darker.

The limiting magnitude of the site (the visual magnitude of the faintest stars detectable with the naked eye at the zenith on a clear night) is estimated from the night sky brightness in magnitudes per square arcsecond from the atlas. Values are converted to a limiting magnitude using equations 17 and 18 from Schaefer 1990 (PASP, 102, 212). The result is only a rough guideline; the actual magnitude of the faintest stars you can see may be higher or lower depending on your age, visual acuity, and observing experience.

User interface notes

Bugs