Astronomy, read from an observer's source
In the places that obeyed Pope Gregory XIII, the day after 4 October 1582 was 15 October. His bull Inter gravissimas of 24 February had removed ten days to bring the spring equinox back to where it stood at the Council of Nicaea in 325, about 21 March. Thursday the 4th was followed by Friday the 15th, and the weekdays did not skip. France followed in December, Britain and its empire not until September 1752, when eleven days went, and Russia in 1918. Anyone who dates a sky must say which calendar the date is in, and astronomers solved the problem by giving up dates. In 1583 the scholar Joseph Scaliger proposed a count of days that began at noon on 1 January 4713 BC in the Julian calendar, and John Herschel brought it into astronomy in his Outlines of Astronomy of 1849.
The astronomy library for DKE, the Deterministic Knowledge Engine, writes some of an observer's work down in DKE Python, and this essay teaches the astronomy by reading it. You record sites, stars, nights and bookings. The library's rules conclude, for each site, the stars that never set, the stars that never rise and the stars it can see, and for each booking the stretch in which the star is up, the night is dark and the instrument is booked. Each answer names the facts it rests on, and correcting a fact withdraws the answers it supported. The second half of the library calculates: dates, sidereal time, altitudes, rising and setting, the Sun, and orbits. Every program below runs as shown once the library's source file has been compiled into a store, and every output is copied from a run. The times are Universal Time, and the angles are in degrees.
1. Dating the sky
The library's jd() turns a calendar date into a Julian date, the day carrying its fraction. It follows the formula of Jean Meeus, the Belgian meteorologist whose Astronomical Algorithms (1991, second edition 1998) is the book many programmers of the sky work from. The formula does what the astronomers do: before the reform it counts in the Julian calendar, and from 15 October 1582 in the Gregorian:
if year < 1582 or (year == 1582 and (month < 10 or (month == 10 and day < 15))):
b = 0So the two dates either side of the gap are a single day apart, and the ten days that never happened are not counted:
import astronomy
print(astronomy.jd(1582, 10, 4.0), astronomy.jd(1582, 10, 15.0))
print(astronomy.jd(2000, 1, 1.5), astronomy.jd(333, 1, 27.5))
print(astronomy.jd(1957, 10, 4.81))2299159.5 2299160.5
2451545.0 1842713.0
2436116.312000 January 1.5 is noon on the first of January 2000, the instant astronomers still call J2000.0, and its Julian date is 2451545.0. 333 January 27.5 is a Julian-calendar date long before the reform, and 1957 October 4.81 is the day, as a fraction, that Sputnik 1 was launched.
2. Sidereal time
A solar day is the time between two noons. In that time the Earth has gone about a degree round its orbit, so it has turned a little more than once relative to the stars, and a day measured by the stars, a sidereal day, is 23 hours 56 minutes 4.09 seconds of ordinary time. A star is back on the meridian then. Sidereal time is the hour angle of the spring equinox, and it says which right ascension is on the meridian now. At Greenwich, Meeus gives it by a polynomial in the Julian date whose main term is the 360.98564736629 degrees the Earth turns in a day of 86,400 seconds:
th = 280.46061837 + 360.98564736629 * (j - 2451545.0) + 0.000387933 * t * t - t * t * t / 38710000gmst() returns it in degrees, lst() adds the longitude of a site, east positive, and hms() writes an angle as the hours, minutes and seconds of the sky's rotation, 15 degrees to the hour. The date is the one Meeus's chapter on sidereal time uses, 1987 April 10 at 19 hours 21 minutes, and Greenwich sidereal time then is 8.5825249 hours, as the IDL Astronomy Library's routine for local sidereal time gives it from the same book:
import astronomy
j = astronomy.jd(1987, 4, 10.0 + 19.35 / 24)
print(j, astronomy.gmst(j))
print(astronomy.hms(astronomy.gmst(j)))
print(astronomy.hms(astronomy.lst(j, 10.75)))2446896.30625 128.73787324433215
[8.0, 34.0, 57.08957863971719]
[9.0, 17.0, 57.089578639717615]Greenwich says 8 hours 34 minutes 57 seconds, and at Oslo, 10.75 degrees east, the sky is 43 minutes further on.
3. Sites and stars
A star is placed by two angles. Its declination is its latitude on the sky, measured from the celestial equator, and its right ascension is its longitude, measured eastward from the spring equinox. A star's magnitude is its brightness on a scale whose history is nearly as long as astronomy's: Ptolemy's Almagest sorted the naked-eye stars into six classes, the brightest first and the faintest sixth, and in 1856 Norman Pogson made the scale exact by letting a first-magnitude star be a hundred times as bright as a sixth. The zero of the modern scale was set from stars of Vega's type, so Vega itself is at 0.03, and Sirius, the brightest star in the night sky, is at -1.46.
The library asks for a site's latitude and longitude, and the faintest magnitude its instrument reaches, and for a star's right ascension, declination and magnitude. The rules then ask one geometric question. At a site of latitude φ the pole stands φ degrees above the horizon, so a star within 90° − φ of it circles without setting. A star the same distance from the other pole never rises. The library computes those two limits for the site, and a rule compares them with each star:
circumpolar_above: real = 90 - latif t.lat >= 0 and s.dec > t.circumpolar_above:
Circumpolar.s.t = TrueHere are five sites and four stars. Polaris is at declination +89.26, 0.74 degrees from the pole; Vega is at +38.78; Sirius at -16.72, and Canopus, the second brightest star in the sky, at -52.70. Oslo is at latitude 59.91 north, Longyearbyen on Svalbard at 78.22, Athens at 37.98, Heraklion on Crete at 35.34 and ESO's Paranal Observatory in the Atacama at 24.63 south:
import astronomy
astronomy.site("oslo", 59.91, 10.75, 6.0)
astronomy.site("longyearbyen", 78.22, 15.65, 6.0)
astronomy.site("athens", 37.98, 23.73, 6.0)
astronomy.site("heraklion", 35.34, 25.14, 6.0)
astronomy.site("paranal", -24.63, -70.40, 6.0)
astronomy.star("polaris", 37.95, 89.26, 1.98)
astronomy.star("vega", 279.23, 38.78, 0.03)
astronomy.star("sirius", 101.29, -16.72, -1.46)
astronomy.star("canopus", 95.99, -52.70, -0.74)
for where in ["oslo", "longyearbyen", "athens", "heraklion", "paranal"]:
print(where, astronomy.circumpolar(where), astronomy.hidden(where))
north = astronomy.Site("oslo").circumpolar_above
south = astronomy.Site("paranal").circumpolar_below
print(north, south)oslo ['polaris', 'vega'] ['canopus']
longyearbyen ['polaris', 'vega'] ['canopus', 'sirius']
athens ['polaris'] ['canopus']
heraklion ['polaris'] []
paranal [] ['polaris']
30.090000000000003 -65.37At Oslo the limit is 30.09 degrees, and Vega, at 38.78, never sets. Sirius is said to be visible from almost everywhere on Earth except north of 73 degrees, and the library gives the same answer where it matters: it is hidden at Longyearbyen, 78 degrees north, and not at Oslo. Canopus is the instructive case. Its declination of -52.70 puts the limit of its rising at 37.3 degrees north, and the encyclopedias give its northern limit in theory as 37 degrees 18 minutes. Athens, at 37.98, is just beyond it and Heraklion, at 35.34, just inside. Refraction lifts a star, and under good conditions it has been seen a little north of the limit; a record from Mount Nemrut in Turkey, at 37 degrees 59 minutes, is the one that is quoted. The library works with the geometric horizon and says what the geometry says.
At Paranal, in the southern hemisphere, the signs are turned round and the limit is -65.37 degrees. No star in the list is within 24.63 degrees of the south celestial pole, so none is circumpolar there, and Polaris, which is nearly at the north pole, never rises.
4. Why a star is up
Polaris is near the pole now, but it has not always been, and the pole does not stay with it. The Earth's axis precesses in a circle that takes about 26,000 years, a motion credited to Hipparchus in the second century BC. Thuban was the pole star around 3000 BC and Kochab from about 1500 BC to AD 500. Polaris will be closest to the pole soon after the year 2100, at about 0.45 degrees. Every coordinate in a catalogue belongs to an epoch, and the declinations above are those of J2000.0.
A conclusion in the library answers why. Ask why Vega is circumpolar at Oslo and the answer is the three facts it rests on: Oslo's latitude, the limit computed from it, and Vega's declination. Then correct Vega's declination to 28.0, below the limit, and the conclusion is withdrawn; correct it back and it returns:
import astronomy
w = why(astronomy.Circumpolar.vega.oslo)
for c in w:
print(str(c.path) + " = " + str(c.value))
astronomy.star("vega", 279.23, 28.0, 0.03)
print(astronomy.circumpolar("oslo"))
astronomy.star("vega", 279.23, 38.78, 0.03)
print(astronomy.circumpolar("oslo"))astronomy.Circumpolar.vega.oslo = True
astronomy.Site.oslo.lat = 59.91
astronomy.Site.oslo.circumpolar_above = 30.090000000000003
astronomy.Star.vega.dec = 38.78
['polaris']
['polaris', 'vega']Nothing was recomputed by hand. The answer followed the fact.
5. Altitude, rising and setting
The altitude of a body above the horizon follows from its declination, the observer's latitude and its hour angle, which is the sidereal time minus its right ascension, so a star on the meridian has an hour angle of 0 and one on the opposite side of the pole 180. semi_arc() solves the same relation the other way round: the hour angle at which the body crosses a given altitude, rising and setting symmetrically about the transit. It answers 0 for a body that never reaches the altitude and 180 for one that never falls below it:
c = (math.sin(math.radians(math.fmod(horizon, 360.0))) - math.sin(p) * math.sin(d)) / (math.cos(p) * math.cos(d))Each angle is taken modulo a full turn first, which math.fmod does exactly, so an angle given past 360 degrees loses no digits on its way to radians.
Polaris stands 60.65 degrees above the horizon at Oslo at its highest and 59.17 at its lowest, and the 59.91 between them is the latitude, which is how a navigator reads latitude from the pole star. rise_set() adds the sidereal time and returns the rising, the transit and the setting of a star on a day, the single transit for one that never sets, and nothing for one that never rises:
import astronomy
print(astronomy.altitude(89.26, 59.91, 0.0), astronomy.altitude(89.26, 59.91, 180.0))
print(astronomy.semi_arc(-16.72, 59.91, 0.0), astronomy.semi_arc(-16.72, 78.22, 0.0), astronomy.semi_arc(38.78, 59.91, 0.0))
t = astronomy.rise_set(101.29, -16.72, 59.91, 10.75, datetime("2026-11-08T12:00:00Z"))
print([str(x)[0:16] for x in t])
print(astronomy.rise_set(101.29, -16.72, 78.22, 15.65, datetime("2026-11-08T12:00:00Z")))
print([str(x)[0:16] for x in astronomy.rise_set(279.23, 38.78, 59.91, 10.75, datetime("2026-11-08T12:00:00Z"))])60.64999999999999 59.17000000000001
58.77387884370289 0.0 180.0
['2026-11-08T22:54', '2026-11-09T02:48', '2026-11-09T06:43']
[]
['2026-11-08T14:42']Sirius is up at Oslo for nearly eight hours, from 22:54 to 06:43 on the night of 8 November, and the library says it never rises at Longyearbyen. The second line has the semi-arcs: 58.8 degrees for Sirius at Oslo, 0 at Longyearbyen and 180 for Vega at Oslo.
6. The Sun and the length of the day
The Sun's declination swings between about -23.4 and +23.4 degrees in a year. sun() gives its right ascension and declination from the low-precision formulas of the Astronomical Almanac, good to about 0.01 degrees from 1950 to 2050, and day_length() uses them to find the Sun's semi-arc to -0.833 degrees: the Sun's upper limb at the horizon after refraction has lifted it. It answers 24 where the Sun never sets and 0 where it never rises:
return 2 * semi_arc(sun(j)[1], latitude, -0.833) / 15l = (280.460 + 0.9856474 * n) % 360At the equator the day is not twelve hours but seven minutes more, because of that 0.833. At Oslo, the summer solstice gives 18 hours 50 minutes and the winter solstice 5 hours 54. The last lines scan the year of 2026 at Tromsø, 69.65 degrees north, for the days the Sun does not set and the days it does not rise:
import astronomy
june = astronomy.jd(2026, 6, 21.5)
december = astronomy.jd(2026, 12, 21.5)
for lat in [0.0, 59.91, 69.65, 78.22]:
print(lat, round(astronomy.day_length(lat, june), 2), round(astronomy.day_length(lat, december), 2))
j0 = astronomy.jd(2026, 1, 1.5)
up = [d for d in range(365) if astronomy.day_length(69.65, j0 + d) == 24.0]
down = [d for d in range(365) if astronomy.day_length(69.65, j0 + d) == 0.0]
print(len(up), str(astronomy.instant(j0 + up[0]))[0:10], str(astronomy.instant(j0 + up[len(up) - 1]))[0:10])
autumn = [d for d in down if d > 200]
spring = [d for d in down if d < 100]
print(len(down), str(astronomy.instant(j0 + autumn[0]))[0:10], str(astronomy.instant(j0 + spring[len(spring) - 1]))[0:10])0.0 12.12 12.12
59.91 18.83 5.9
69.65 24.0 0.0
78.22 24.0 0.0
69 2026-05-18 2026-07-25
48 2026-11-28 2026-01-14The library has the midnight sun at Tromsø for 69 days, from 18 May to 25 July. Wikipedia gives the Sun above the horizon from about 18 May to 26 July, which differs by a day at the end, and from 19 May to 27 July in another place in the same article. The page adds that the mountains round the city hide the Sun for a few days at each end, so that it is seen from about 21 May to 21 July. The polar night is 48 days from 28 November to 14 January, which is the figure and the dates the article gives for the continuous night. At Longyearbyen, 78 degrees, the Sun is above the horizon from 18 April to 24 August, the article says, and below it from 27 October to 15 February. It does not appear in the town until about 8 March, because of the mountains, and the library, which has no mountains, does not know that.
7. Nights and bookings
An observer needs a night as well as a sky. Astronomers divide the twilight by how far the Sun is below the horizon: civil twilight to 6 degrees, nautical to 12, and astronomical to 18, beyond which it is night. The library does not compute nights; you write them, and it works with what you write. This program computes the dark of 8 November at Oslo with the library's own scripts, taking the hour angle at which the Sun is 18 degrees below the horizon from semi_arc() and placing it round the Sun's transit, and writes it as a night:
for b in Booking, n in Night, u in Up:
if n.site == b.site and u.site == b.site and u.star == b.star:
b.observable = TrueA booking is an instrument at a site for a star over a stretch of time. up() writes the stretch a star is above the site's horizon, and a rule concludes that the booking is observable for the overlap of the booking, the night and the star's time up, and for no more. window() reports the overlap as its first and last instants. Vega is up all day at Oslo, so its booking from 15:00 to 22:00 is cut only by the dusk. Sirius's, from 18:00 to 05:00, is cut at both ends, by its rising at 22:54 and by the dawn. Canopus never rises at Oslo, so its booking is never observable:
import astronomy
day = datetime("2026-11-08T00:00:00Z")
s = astronomy.sun(astronomy.jd(2026, 11, 8.5))
noon = astronomy.rise_set(s[0], s[1], 59.91, 10.75, day)[1]
w = duration("P1D") * (astronomy.semi_arc(s[1], 59.91, -18.0) / 360.0)
dusk = noon + w
dawn = noon + duration("P1D") - w
print(str(dusk)[0:16], str(dawn)[0:16])
astronomy.night("n1", "oslo", dusk, dawn)
astronomy.up("vega", "oslo", dusk)
astronomy.book("b1", "oslo", "vega", datetime("2026-11-08T15:00:00Z"), datetime("2026-11-08T22:00:00Z"))
print([str(x)[0:16] for x in astronomy.window("b1")])
astronomy.up("sirius", "oslo", dusk)
astronomy.book("b2", "oslo", "sirius", datetime("2026-11-08T18:00:00Z"), datetime("2026-11-09T05:00:00Z"))
print([str(x)[0:16] for x in astronomy.window("b2")])
astronomy.book("b3", "oslo", "canopus", datetime("2026-11-08T18:00:00Z"), datetime("2026-11-09T05:00:00Z"))
print(astronomy.window("b3"))
print(astronomy.observable_at(datetime("2026-11-08T20:00:00Z")), astronomy.observable_at(datetime("2026-11-09T01:00:00Z")), astronomy.observable_at(datetime("2026-11-09T04:45:00Z")))2026-11-08T17:30 2026-11-09T04:31
['2026-11-08T17:30', '2026-11-08T22:00']
['2026-11-08T22:54', '2026-11-09T04:31']
[]
['b1'] ['b2'] []The first booking is observable from the end of twilight, at 17:30, to its own end at 22:00. The second begins at 22:54, when Sirius rises, and ends at 04:31 with the night, though the booking runs to 05:00. The third is empty. The last line reads at an instant, which is what a night's duty astronomer wants to know: at 20:00 the first booking is observable, at 01:00 the second, and at 04:45 neither.
Observatories do this on a larger scale and for a living. The Very Large Telescope, four 8.2-metre telescopes, stands at Paranal, at 24 degrees 37 minutes south and 2,635 metres, in an observatory the European Southern Observatory established in 1999, and which stars it can reach on a given night is the question of this section.
8. When two sources disagree
Two catalogues will disagree: they are written for different epochs, or by different hands. A second source writes Paranal's latitude, 24 degrees 37 minutes 38 seconds south, as -24.6272 where the first rounded it to -24.63. A second catalogue, invented for the demonstration, puts Canopus at -20.0 degrees. The library does not choose between the two. It reasons about both, no query refuses, and unsettled() names what they disagree about:
import astronomy
update(astronomy.Site("paranal").lat, -24.6272, "second source")
update(astronomy.Star("canopus").dec, -20.0, "second catalogue")
print(astronomy.unsettled())
print(astronomy.hidden("oslo"), astronomy.visible("oslo"))['site paranal lat', 'star canopus dec']
['canopus'] ['canopus', 'polaris', 'sirius', 'vega']Canopus is hidden at Oslo by one declination and visible by the other, so it appears in both lists, and the fact that is in dispute is named. Only up(), which needs one declination to compute a passage, would be refused until the dispute is settled.
9. Kepler's equation
A planet in a circular orbit moves at a constant rate. In an ellipse it does not, and the question of where it is at a given time is Kepler's, in his Astronomia nova of 1609, where he found that Mars moves on an ellipse and sweeps equal areas in equal times. If the mean anomaly M is the fraction of the period elapsed, written as an angle, the eccentric anomaly E, the angle that fixes the body's place on the circle round its ellipse, satisfies E − e·sin E = M, where e is the eccentricity. The equation cannot be solved for E in closed form, and the usual way today is Newton's method. kepler() starts from M, or from π when e is above 0.8 and a start at M might stall, and takes steps until one changes E by less than 10⁻¹⁵:
step = (x - e * math.sin(x) - m) / (1 - e * math.cos(x))if e > 0.8:
x = math.pikepler() returns E and the number of steps it took, and true_anomaly() turns E into the angle seen from the focus. The Earth's orbit has an eccentricity of about 0.0167, and Halley's comet 0.967:
import astronomy
print(astronomy.kepler(1.0, 0.0167))
print(astronomy.kepler(0.1, 0.967))
k = astronomy.kepler(1.0, 0.5)
print(k, astronomy.true_anomaly(k[0], 0.5))[1.0141790871647136, 4.0]
[0.7802277443640641, 8.0]
[1.4987011335178484, 6.0] 2.030806214849156Four steps for the Earth and eight for the comet, which spends most of its time far from the Sun. Edmond Halley argued in his Synopsis of the Astronomy of Comets of 1705 that the comets of 1531, 1607 and 1682 were one body returning about every 76 years, and predicted its return in 1758. He died in 1742; it was first seen again on Christmas Day 1758 by Johann Georg Palitzsch, a German farmer and amateur astronomer. It passed the Sun in 1986, and its orbital period is about 74.7 years, so it is due again on 28 July 2061.
10. Two bodies, integrated
Kepler's equation gives the place on the ellipse at once. A computer can also take the orbit a step at a time, from the force, and that is the only way when other bodies pull. two_body() does it for one body about a mass in units where the mass times the gravitational constant is 1 and the semi-major axis is 1, starting at perihelion, with the classical fourth-order method that Carl Runge and Wilhelm Kutta worked out about 1900, which takes four slopes in a step and weights them:
x = x + dt / 6 * (vx + 2 * u2 + 2 * u3 + u4)It returns the energy at the start, the energy at the end, and the distance between where the steps put the body and where Kepler's equation puts it. The energy of a two-body orbit does not change, so how far it drifts is an error that needs no answer book. A near-circular orbit at 200 steps to the orbit comes through almost exactly; Halley's comet does not:
import astronomy
print(astronomy.two_body(0.0167, 1.0, 200))
print(astronomy.two_body(0.967, 1.0, 200))
print(astronomy.two_body(0.967, 1.0, 2000))
print(astronomy.two_body(0.967, 1.0, 10000))[-0.5, -0.5000000026954337, 1.6154076665234124e-07]
[-0.5, 31.29799145995318, 49.82083970051012]
[-0.5, -0.6157642907074433, 1.237684286453551]
[-0.5, -0.5000297785903172, 0.006286905305969945]The first line is the Earth's orbit, with an energy that has moved in the ninth decimal and a position wrong by under a millionth of the radius. The second is the comet at 200 steps to the orbit, whose energy has gone from -0.5 to +31.3, an unbound orbit, a comet thrown out of the solar system by a numerical error. At 2,000 steps the orbit is bound again but badly wrong. It takes 10,000 steps an orbit to bring the position within a hundredth. The step is fixed, and a comet asks for small steps near the Sun and large ones far from it. Programmers of orbits for that reason let the step change; this library takes the same step everywhere.
11. What the library leaves out
Rising and setting are for the geometric horizon unless a horizon is given; day_length() is the only script that includes refraction, which lifts a body by about half a degree. A star's coordinates are the ones you write, for the epoch you choose, and nothing in the library applies precession, nutation, aberration, proper motion or parallax. A site is a point with no horizon of hills and no altitude. Visible asks only whether the star is ever above the horizon and bright enough, and the library says nothing of the Moon, the weather, or the time of year. The Sun is the low-precision Sun, there is no Moon or planet, and the orbits are one body about a fixed mass. A night is whatever you write. These are the textbook cases, and each has an answer known some other way.
Closing
The sky was dated with a count of days because calendars could not agree on the days, and its oldest measurements, the six classes of stars and the circle of precession, were the work of people who wrote down what they saw, and later of people who could compute what they would see. The library is on that side of the work. It counts the days without the calendar, says which stars a site can ever see and which it can never, lays a booking against a night and a rising and reports the stretch that survives, and when two sources disagree it says so. And because each conclusion names what it rests on, a reader who doubts one can ask why, and follow the answer down to the latitude, the declination or the night it came from.