The starting and finishing points of astronomical calculations are usually angles, expressed as degrees, minutes and seconds of arc, or times, measured in hours, minutes, and seconds. Yet computers often cannot handle these quantities as they stand. They must first be converted to decimal degrees or hours before the computation can begin, and then the result converted back into the more-familiar minutes and seconds form. Consider, for example, the problem of finding the time corresponding to 1 hour 37 minutes before 3 : 25 in the afternoon. We first convert the times to their decimal forms, 1.616667 hours and 3.416667 hours, then subtract the former from the latter to get 1.800 hours, and finally convert back to minutes and seconds: 1 : 48 pm.

Your computer can deal with the problem of converting angles or times between decimal and minutes/seconds forms with the aid of the subroutine MINSEC. The direction in which the conversion is carried out is controlled by the switch SW(1). When SW(1) has the value +1, then the quantity represented by the input argument X (decimal degrees/hours) is converted into degrees/hours, minutes, and seconds which are returned by the output arguments XD, XM, XS respectively. XD and XM represent whole numbers while XS represents the seconds correct to two decimal places. The sign of X is returned by the string variable S$ which has the value ‘ + ’ unless X is negative when S$ becomes ‘―’. The real variable SN also carries this information since SN = +1 if X is positive, SN = 0 if X is zero, and SN = −1 if X is negative.

The prediction of the position of a planet, comet, or any other member of the Solar System, depends upon knowing its orbital elements. The orbit is completely defined by these elements (in the absence of perturbations) and hence the position of the body can be calculated for any moment in the future or the past. Sometimes, however, you may not know the elements, particularly if you have been fortunate enough to have discovered a new comet. In that case, you can measure the position of the comet on at least three separate occasions as far apart as possible and then apply one of a number of methods to deduce the set of elements which are consistent with your observations. Unfortunately, the standard methods are neither particularly easy to apply, nor to understand, so I have adopted a more obvious, if perhaps rather flat-footed, approach to the problem which makes use of the mindless number-crunching ability of a computer.

The method is best understood by first considering a simpler problem in one dimension only. Suppose that we have a mathematical relationship between the variables VO and P such that VO is a function of P i.e. given any value of P we can immediately calculate the value of VO.

The path of any of the major planets in its elliptical orbit around the Sun may be defined in terms of six orbital elements:the mean longitude, the longitude of the perihelion, the eccentricity, the inclination, the longitude of the ascending node, and the length of the semi-major axis. Of these, only the last is a true constant. The other elements, being referred to the planes of the ecliptic and the equator, reflect the slow changes caused by the influences of other members of the Solar System. It is therefore necessary to calculate the elements for a given date, and the orbit so defined will then be the correct mean orbit of the planet at that time. The true orbit may be further perturbed by short- and longperiod terms which must be considered if high accuracy is required (see subroutine PLANS, 4500).

The orbital elements are expressed by polynomials of the form

A0 + A1T + A2T2 + A3T3

where A0, Al, etc. are constants and T is the number of Julian centuries of 36525 Julian days since 1900 January 0.5. Routine PELMENT calculates the values of the orbital elements for all the planets on a given date (DY, MN, YR) and returns them in the array PL (8,9). This array must be declared before calling the routine. PELMENT also returns the mean daily motion in longitude, the angular diameter at a distance of 1 AU from the Earth, and the standard visual magnitude.

The combined gravitational fields of the Sun and the Moon acting on the nonspherical Earth cause the direction of the Earth's rotation axis to gyrate slowly with a period of about 25 800 years. This effect, called precession, is calculated by routines PRCESS1 (2500), and PRCESS2 (2600). Superimposed on the regular motion there are also small additional periodic terms caused by the varying distances and relative directions of the Moon and Sun, which continuously alter the strength and direction of the gravitational field. This slight wobbling motion is called nutation, and it must be taken into account when accuracies of better than half an arcminute are required. The routine given here calculates the effects of nutation on the ecliptic longitude (DP) and on the obliquity of the ecliptic (DO). It has an intrinsic accuracy (reduced by lack of precision in your machine) of about 1 arcsecond.

The calculations are made by the routine for the date input as DY (days), MN (months), and YR (years). Execution is controlled by the flag FL(6). If this flag is set to 0 before calling NUTAT, execution of the routine continues normally. It is set to 1 at line 1865. Subsequent calls to NUTAT with FL(6) = 1 do not result in any new calculations, but control is returned immediately to the calling program.

Routine RISET (3100) calculates the times and azimuths of rising and setting of a celestial body given its right ascension and declination as input parameters. We may therefore combine it with the routines MOON (6000) and EQECL (2000) to find the circumstances of moonrise and moonset at any place on the Earth's surface for any given calendar date. There is a snag, however, as the Moon's position is constantly changing by about half a degree of ecliptic longitude per hour and, unless we already know the times of rising and setting, we cannot find the correct position of the Moon in advance. If we use the position calculated for a set time, say midnight, we may be as much as one hour out in the calculation of moonrise and moonset. We must therefore adopt an iterative procedure, refining an initial crude estimate to achieve the required accuracy, just as in SUNRS (3600).

Routine MOONRS first calculates the right ascension and declination of the Moon at local midday of the given date, making allowance for nutation (NUTAT, 1800). Then the universal times (UTs) and azimuths of the corresponding moonrise and moonset are found via RISET (3100) and TIME (1300) with the corrections for parallax, refraction and the angular size of the Moon's disk incorporated in the parameter DI. The Moon's position is now recalculated at each of these approximations to moonrise or moonset and the whole procedure repeated twice more to find times of rising and setting accurate to within one minute, and azimuths to within about 1 arcminute.

The starting point for calculations involving a body in an elliptical orbit is often the mean anomaly, AM. This is the angle moved by a fictitious body in a circular orbit of the same period as the real body, the angle being reckoned in the same sense as the direction of motion of the real body from the point of closest approach (the periapsis).The quantity needed is the true anomaly, AT, which measures the angle moved by the real body since periapsis, and it is related to AM through the eccentric anomaly, AE, by Kepler's equation.

AE - (EC × SIN(AE)) = AM,

where EC is the eccentricity of the orbit. Unfortunately, this equation is not easily solved, but the solution can be approximated by a trigonometric expansion called the equation of the centre.If EC is less than about 0.1 and high precision is not required, the first term of the expansion may suffice, giving

AT = AM + (2 × EC × SIN(AM)),

where AT and AM are expressed in radians.

For more accurate work, the equation must be solved explicitly for AE, and then AT calculated from

TAN(AT/2) = ((1 + EC)/(1 - EC))0.5 × TAN(AE/2).

The routine given here solves Kepler's equation by an iterative method in which an approximate solution for AE is repeatedly refined until the error between (AE - (EC x SIN(AE))) and AM is less than a given error (10-6 radians).

The eclipse pictures printed in Figure 7, which show lunar and solar eclipses at various stages, were produced by subroutine DISPLAY listed later in this section. Unfortunately, it is not possible to write code which will work equally well on every computer since the mechanics of producing graphics displays varies so much from machine to machine. You will certainly have to modify DISPLAY to suit the particular characteristics of your own computer. To help you do so, I have collected all the graphics commands into a number of function calls; for example, FNCLGS clears the graphics screen to black. You can use DISPLAY just as it is provided you can define the functions so that they work on your machine. Alternatively, you may replace the functions with calls to procedures or subroutines to do the same thing. Thus, if your BASIC includes the graphics command ‘HOME’, for example, which clears the screen, you could substitute ‘GOSUB 9000’ for ‘X = FNCLGS’ and include the extra subroutine

9000 Home : RETURN

at the end. Obviously, more complex procedures will require several lines of code each.

The graphics screen on my computer measures 392 pixels in the horizontal, or X direction, and 256 in the vertical, or Y direction. You can think of the screen as a large sheet of graph paper, with the point X = 0, Y = 0 at the bottom left-hand corner of the screen.

The right ascension and declination of a celestial body change slowly with time because of luni-solar precession, a gyrating motion of the Earth's axis caused by the gravitational effects of the Sun and the Moon on the equatorial bulge of the Earth. The Earth's axis describes a cone about a line perpendicular to the ecliptic through the centre of the Earth with half-angle approximately 23.5 degrees, completing one circuit in about 25 800 years. The change is thus about 50 arcseconds per year and this must be allowed for when comparing equatorial coordinates at two different epochs.

PRCESS1

Two subroutines are given here which will correct equatorial coordinates for the effects of precession between one epoch and another. The first of these, PRCESS1 (2500), uses an approximate method which is simple and compact, but which nevertheless gives good results. You should find it sufficient for most purposes. The corrections, XI and Yl, which have to be made to the right ascension, X, and declination, Y, measured at epoch A, to find the corresponding values at epoch B, are well approximated by the expressions:

X1 = (MP + (NP × SIN(X) × TAN(Y))/15) × NY

Y1 = NP × COS(X) × NY,

where NY is the number of years between the two epochs, and MP and NP are ‘constants’ which change slowly with time.

The basis of all civil timekeeping on the Earth is Universal time, UT. It is a measure of time which conforms, within a close approximation, to the mean daily motion of the Sun. Roughly speaking, the time is 12 : 00 noon UT when the Sun crosses the meridian at Greenwich, and 24 hours of UT is the time between two successive passages of the Sun across that meridian. However, the Sun is not really a very good timekeeper. Although the Earth spins about its polar axis at a nearly-constant rate (there are small irregular perturbations), the position of the Sun in the sky changes daily due to the orbit of the Earth about the Sun. If that orbit were perfectly circular, and if the Earth's axis were exactly perpendicular to the plane of the orbit (the ecliptic), then the daily passage of the Sun would be regular throughout the year. In reality, the Earth's orbit is slightly elliptical so that the speed of the Earth varies throughout the year, and the Earth's axis is tilted by about 23.5 degrees from the perpendicular to the plane of the ecliptic. Both these effects cause the time measured by a sundial, true solar time, to run fast during some parts of the year and slow during others by up to 16 minutes with respect to a uniform clock. The difference between the solar time and the uniform clock is called the ‘equation of time’, and you would need to know its value if you wished to tell the time accurately from a sundial.

The equatorial coordinates of members of the Solar System calculated in later subroutines are those appropriate to an observer situated at the centre of the Earth (geocentriccoordinates). When the celestial body in question is at a very large distance away, the same coordinates apply to observations made from the Earth's surface. However, objects within the Solar System are sufficiently close that their apparent positions (topocentriccoordinates) change slightly depending on the exact vantage point. The Sun, for example, may appear displaced by as much as 8 arcseconds from its geocentric position, and the Moon by as much as 1 degree.

Routine PARALLX calculates the apparent hour angle (P) and declination (Q) for an observer at a given geographical latitude (GP) and height above sealevel (HT), when the geocentric hour angle (X), declination (Y), and equatorial horizontal parallax (HP) are known. It will also perform the reverse operation, calculating the geocentric coordinates from the apparent position. The correction for parallax to the geocentric hour angle (DX) is the same as that to the geocentric right ascension. Having run the routine, the apparent right ascension, RA1, may be found from

RA1 = RA - DX,

where RA is the geocentric right ascension.

Execution of PARALLX is controlled by the flag FL(12). This should be set to 0 before calling the routine for the first time, or whenever new values of GP or HT are specified.

The subroutines described earlier in this book for converting between the various coordinate systems are quite satisfactory for normal use where you have a requirement for one particular conversion. However, it is sometimes convenient to have a more general routine to hand which can be used to convert from one system to any other system, at the flick of a (software) switch. We can use matrices for this purpose, ordered sets of numbers set out in rows and columns. The manipulation of the matrices is always the same; to change from one system to another you merely have to change the numbers in the matrices.

Subroutine GENCON is able to produce and manipulate just four different matrices, appropriate for conversions between azimuth/altitude and hour angle/declination, hour angle/declination and right ascension/declination, right ascension/declination and ecliptic longitude/latitude, and right ascension/ declination and galactic longitude/latitude. However, you can use it to convert coordinates from any of the systems to any other system by serial conversion. For example, suppose that you wish to know the azimuth and altitude of a pulsar whose galactic coordinates you have to hand. You would first convert to right ascension/declination, then to hour angle/declination, and finally to azimuth/altitude. GENCON will do this triple conversion automatically.

An eclipse occurs whenever the Moon's apparent position coincides with that of the Sun (solar eclipse) or the Earth's shadow (lunar eclipse). In the former case, the Moon is seen to pass in front of the Sun and obscure it. In the latter case, the bright Moon fades as it passes into the shadow of the Earth, until it is dimly illuminated only by the reddened light refracted through the Earth's atmosphere if the eclipse is total. A solar eclipse can happen only at new moon and a lunar eclipse at full moon. There is not an eclipse on every such occasion, however, since the Moon's orbit is inclined at an angle of about 5 degrees to the plane of the ecliptic. This means that the Moon must also be within about 18.5 degrees of one of its nodes, limiting the total number of eclipses possible each year to a maximum of seven.

Subroutine ECLIPSE (7000) calculates all the circumstances of a lunar or solar eclipse occurring near a given date. This date is input as usual via DY, MN, YR, and the type of eclipse via the string variable ET$, set to 'L’ for a lunar eclipse and to ‘S’ for a solar eclipse. The observer's position on the Earth is also input via the variables GL (geographical longitude), GP (geographical latitude) and HT (height above sea-level).

Comets are members of our Solar System, usually with highly elongated orbits, which become visible near the Sun. They have bright heads and diffuse tails of variable length which always point away from the Sun. Comets can be divided into two categories: those which are gravitationally bound to the Sun like the planets and travel in elliptical orbits, the periodiccomets, and those which do not seem to be bound to the Sun but appear once and then shoot off into space never to be seen again, the paraboliccomets. Routine PCOMET deals with the latter sort; you can use routine ELOSC (7500) for periodic comets.

The approximate position of a parabolic comet can be found from its parabolic orbital elements. As the name suggests, the orbit is a parabola and needs one less element to define it than does an ellipse. Subroutine PCOMET (7900) performs this task, taking the five parabolic orbital elements and calculating the ecliptic coordinates of the comet for a given instant. This is input in the usual way via DY, MN, YR, returning with the ecliptic longitude, EP</b., the latitude, BP, and the distance of the comet from the Earth, RH. These quantities have all been corrected for light-travel time. The routine also returns the instantaneous values (i.e. not corrected for light-travel time) of the heliocentric ecliptic longitude, L0, latitude, SO, radius vector, P0, and distance from Earth, V0.

The problem of converting between the hour angle and the right ascension of a celestial object at a given local sidereal time is so straightforward as to make it hardly seem worthwhile devoting a separate subroutine to it. Nevertheless, it crops up sufficiently often that one quickly becomes tired of writing the few lines of code needed every time into a new program, and so I have included the routing HRANG here. The equation is the same for conversion in either direction, there being no need for a switch to specify it. The right ascension or hour angle (hours) is input via the variable X and the corresponding hour angle or right ascension (hours) returned by P. The geographical longitude, GL, must be specified in degrees (W negative) and the Greenwich sidereal time, SG, in hours. The routine also returns the local sidereal time, LS, in hours.

Formulae

1. P = [LS − X]24

2. LS = [SG + (GL/15)]24

3. P = hour angle or right ascension

given

1. X = right ascension or hour angle

2. SG = Greenwich sidereal time

3. GL= geographical longitude

[]24 indicates reduction to the range 0-24 by the addition or subtraction of multiples of 24.

Details of HRANG

Called by GOSUB 1600.

Converts X, the hour angle or right ascension, into P, the corresponding right ascension or hour angle, at SG, the Greenwich sidereal time, all quantities in hours.The geographical longitude, GL, must also be specified in degrees(W negative, E positive). The local sidereal time in hours is returned by LS.

Many problems in computational astronomy take a date and a time as their starting points. For example, we may wish to calculate the position of Venus at 4 : 30 pm on 17th August 1938, or the phase of the Moon at midnight on the 19th May 1997. These dates and times are reckoned according to some nationally-agreed calendar which ascribes numbers to each instant so that it may be identified uniquely. If we say that a particular event occurred at 11 o'clock in the morning on 17th August 1938, we usually mean that it occurred at 11 hours after the beginning of the 17th day after the beginning of the 8th month after the beginning of the 1938th year after the adopted instant of the birth of Christ. At least, this is what we mean if we are using the Gregorian calendar, which is usually the case.

In order that we may make use of the given date in a computer program, it must first be converted to a single number which can be handled easily by the machine. We must choose some particular instant as our starting point, and then count the days logically from that moment. There are several obvious choices. We could, for example, choose the agreed instant of the birth of Christ, as the Gregorian calendar seems to do.