observability.date is a tool for planning astronomical observations. It 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. To calculate the parameters I want the website to provide, I relied heavily on Astronomical Algorithms by Jean Meeus. I also referred to a few other websites and resources.
I have tested observability.date in several time zones, but only on Linux machines running Google Chrome and Mozilla Firefox. 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.
This is the main functionality of observability.date, and also just about the simplest part to code. Calculating the altitude and azimuth of a celestial object from its Right Ascension and Declination just requires a calculation of the local sidereal time, 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, neglecting the correction for nutation which amounts to only 1.1s at most.
In calculating the altitude of objects other than the Sun and the Moon, I neglect the effect of refraction, which delays their setting and advances their rising by a couple of minutes. As the tool is for finding the best times to observe things, their movements very close to the horizon don't need to be treated carefully.
When an object's coordinates are resolved from its name, observability.date 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 is a small effect but still amounts to 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.
From the altitude of the object I calculate the airmass that it is at. The simplest equation for airmass assumes a plane parallel atmosphere and is X=sec(z), where z is the zenith distance. I correct for the curvature of the Earth using the equation of Rozenberg (1966), which is
X = (cos(z) + 0.25e-11cos(z))-1
To calculate sunset and sunrise times, observability.date 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 observability.date 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 observability.date 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.
This is by far the most complicated thing that observability.date calculates. 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 somewhat simplified set of calculations in chapter 45 of Astronomical Algorithms, with 59 periodic terms in longitude and 30 in latitude. I simplify this further by taking the first 13 periodic terms in the Moon's longitude and the first 10 in latitude.
An exact calculation of moonrise and moonset needs to be done iteratively, as its RA and dec are changing constantly. observability.date 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 couple of 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.
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. observability.date calculates the parallactic angle using the following equation, and displays it in the dynamic display on the daily chart.
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.
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.
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.
I've compared the numbers from observability.date 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 observability.date (calculated on 16 Jan 2018, internal version ID a38519), with results from the ING's staralt, Heavens Above, and ESO's skycalc.
| observability.date | heavens above | eso skycalc | staralt | |
| Sunset | 18:15:34 | 18:08 | 18:16 | 18:15 |
| End of civil twilight | 18:31:45 | 18:32 | ||
| End of nautical twilight | 18:59:58 | 19:01 | 19:01 | |
| End of astronomical twilight | 19:27:45 | 19:28 | 19:28 | 19:27 |
| Start of astronomical twilight | 06:05:20 | 06:06 | 06:06 | 06:05 |
| Start of nautical twilight | 06:33:06 | 06:33 | 06:33 | |
| Start of civil twilight | 07:01:17 | 07:01 | ||
| Sunrise | 07:17:27 | 07:26 | 07:18 | 07:18 |
| Length of night (sunset-sunrise) | 13:02 | 13:18 | 13.0 | 13:03 |
| Length of astronomical night | 10:38 | 10:38 | 10.6 | 10:38 |
| LST at civil midnight | 18:30 | 18:30 | ||
| Moon altitude at midnight | -58.9° | |||
| Distance from object | 124.5° | 123 | ||
| Illuminated fraction | 0.13 | 0.14 | 0.122 | 0.12 |
| Days past full | 12 | |||
| Sets | 15:12:15 | 15:01 | 15:09 | |
| Rises | 04:22:46 | 04:36 | 04:28 | 04:40* |
| object max alt | 48.6 | 48* | ||
| airmass | 1.33 | 1.35* | ||
| parallactic angle at midnight | 86 | 79 | ||
| hours above 30 degrees in dark: | 4h18m* | |||
| * = 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. There are three or four likely causes of this. Firstly, the precision of the lunar calculations; observability.date uses only the first 13 periodic terms in lunar longitude tabulated in Meeus. The details of the calculations for the other programs are not documented, as far as I am aware. Secondly, the lunar position may be geocentric rather than topocentric; lack of available documentation prevents a conclusion from being drawn on that. Thirdly, altitude may not be corrected for; and finally the elevation of the observatory may not be corrected for. This is certainly the case for Heavens Above, which explicitly states that it uses a horizon altitude of -0.8 degrees, which places the refracted upper edge of the Sun or Moon onto the horizon at sea level.
The heliocentric velocity correction calculated by observability.date 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; I've gone with the latter).
Most of the equations used in observability.date 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. observability.date 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: observability.date 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.
Information is returned about the darkness of the night sky at the selected observing location. This information is taken from the World Atlas of the Night Sky Brightness, as updated in 2006. 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.
d toggles between normal and dark appearance