4 EXPLANATION AND EXAMPLES

To guide the writer of positional-astronomy applications software, this final chapter puts the SLALIB routines into the context of astronomical phenomena and techniques, and presents a few “cookbook” examples of the SLALIB calls in action. The astronomical content of the chapter is not, of course, intended to be a substitute for specialist text-books on positional astronomy, but may help bridge the gap between such books and the SLALIB routines. For further reading, the following cover a wide range of material and styles:

• Explanatory Supplement to the Astronomical Almanac, ed. P. Kenneth Seidelmann (1992), University Science Books.
• Vectorial Astrometry, C. A. Murray (1983), Adam Hilger.
• Spherical Astronomy, Robin M. Green (1985), Cambridge University Press.
• Spacecraft Attitude Determination and Control, ed. James R. Wertz (1986), Reidel.
• Practical Astronomy with your Calculator, Peter Duffett-Smith (1981), Cambridge University Press.

Also of considerable value, though out of date in places, are:

• Explanatory Supplement to the Astronomical Ephemeris and the American Ephemeris and Nautical Almanac, RGO/USNO (1974), HMSO.
• Textbook on Spherical Astronomy, W. M. Smart (1977), Cambridge University Press.

Only brief details of individual SLALIB routines are given here, and readers will find it useful to refer to the subprogram specifications elsewhere in this document. The source code for the SLALIB routines (available in both Fortran and C) is also intended to be used as documentation.

4.1 Spherical Trigonometry

Celestial phenomena occur at such vast distances from the observer that for most practical purposes there is no need to work in 3D; only the direction of a source matters, not how far away it is. Things can therefore be viewed as if they were happening on the inside of sphere with the observer at the centre – the celestial sphere. Problems involving positions and orientations in the sky can then be solved by using the formulae of spherical trigonometry, which apply to spherical triangles, the sides of which are great circles.

Positions on the celestial sphere may be specified by using a spherical polar coordinate system, defined in terms of some fundamental plane and a line in that plane chosen to represent zero longitude. Mathematicians usually work with the co-latitude, with zero at the principal pole, whereas most astronomical coordinate systems use latitude, reckoned plus and minus from the equator. Astronomical coordinate systems may be either right-handed (e.g. right ascension and declination $\left[\alpha ,\delta \right]$, Galactic longitude and latitude $\left[lII,bII\right]$) or left-handed (e.g. hour angle and declination $\left[h,\delta \right]$). In some cases different conventions have been used in the past, a fruitful source of mistakes. Azimuth and geographical longitude are examples; azimuth is now generally reckoned north through east (making a left-handed system); geographical longitude is now usually taken to increase eastwards (a right-handed system) but astronomers used to employ a west-positive convention. In reports and program comments it is wise to spell out what convention is being used, if there is any possibility of confusion.

When applying spherical trigonometry formulae, attention must be paid to rounding errors (for example it is a bad idea to find a small angle through its cosine) and to the possibility of problems close to poles. Also, if a formulation relies on inspection to establish the quadrant of the result, it is an indication that a vector-related method might be preferable.

As well as providing many routines which work in terms of specific spherical coordinates such as $\left[\alpha ,\delta \right]$, SLALIB provides two routines which operate directly on generic spherical coordinates: sla_SEP computes the separation between two points (the distance along a great circle) and sla_BEAR computes the bearing (or position angle) of one point seen from the other. The routines sla_DSEP and sla_DBEAR are double precision equivalents. As a simple demonstration of SLALIB, we will use these facilities to estimate the distance from London to Sydney and the initial compass heading:

(The result is 17011 km, $61^{\circ }$.)

The routines sla_SEPV, sla_DSEPV, sla_PAV, sla_DPAV are equivalents of sla_SEP, sla_DSEP, sla_BEAR and sla_DBEAR but starting from vectors instead of spherical coordinates.

4.1.1 Formatting angles

SLALIB has routines for decoding decimal numbers from character form and for converting angles to and from sexagesimal form (hours, minutes, seconds or degrees, arcminutes, arcseconds). These apparently straightforward operations contain hidden traps which the SLALIB routines avoid.

There are five routines for decoding numbers from a character string, such as might be entered using a keyboard. They all work in the same style, and successive calls can work their way along a single string decoding a sequence of numbers of assorted types. Number fields can be separated by spaces or commas, and can be defaulted to previous values or to preset defaults.

Three of the routines decode single numbers: sla_INTIN (integer), sla_FLOTIN (single precision floating point) and sla_DFLTIN (double precision). A minus sign can be detected even when the number is zero; this avoids the frequently-encountered “minus zero” bug, where declinations etc. in the range $0^{\circ }$ to $-1^{\circ }$ mysteriously migrate to the range $0^{\circ }$ to $+1^{\circ }$. Here is an example (in Fortran) where we wish to read two numbers, an integer IX and a real, Y, with IX defaulting to zero and Y defaulting to IX:

Two additional routines decode a 3-field sexagesimal number: sla_AFIN (degrees, arcminutes, arcseconds to single precision radians) and sla_DAFIN (the same but double precision). They also work using other units such as hours etc. if you multiply the result by the appropriate factor. An example Fortran program which uses sla_DAFIN was given earlier, in section 1.2.

SLALIB provides four routines for expressing an angle in radians in a preferred range. The function sla_RANGE expresses an angle in the range $\pm \pi$; sla_RANORM expresses an angle in the range $0-2\pi$. The functions sla_DRANGE and sla_DRANRM are double precision versions.

Several routines (sla_CTF2D, sla_CR2AF etc.) are provided to convert angles to and from sexagesimal form (hours, minute, seconds or degrees, arcminutes and arcseconds). They avoid the common “converting from integer to real at the wrong time” bug, which produces angles like $24h59m$$59$$s$$.999$. Here is a program which displays an hour angle stored in radians:

4.2 Vectors and Matrices

As an alternative to employing a spherical polar coordinate system, the direction of an object can be defined in terms of the sum of any three vectors as long as they are different and not coplanar. In practice, three vectors at right angles are usually chosen, forming a system of Cartesian coordinates. The x- and y-axes lie in the fundamental plane (e.g. the equator in the case of $\left[\alpha ,\delta \right]$), with the x-axis pointing to zero longitude. The z-axis is normal to the fundamental plane and points towards positive latitudes. The y-axis can lie in either of the two possible directions, depending on whether the coordinate system is right-handed or left-handed. The three axes are sometimes called a triad. For most applications involving arbitrarily distant objects such as stars, the vector which defines the direction concerned is constrained to have unit length. The x-, y- and z-components can be regarded as the scalar (dot) product of this vector onto the three axes of the triad in turn. Because the vector is a unit vector, each of the three dot-products is simply the cosine of the angle between the unit vector and the axis concerned, and the x-, y- and z-components are sometimes called direction cosines.

For some applications, including those involving objects within the Solar System, unit vectors are inappropriate, and it is necessary to use vectors scaled in length-units such as AU, km etc. In these cases the origin of the coordinate system may not be the observer, but instead might be the Sun, the Solar-System barycentre, the centre of the Earth etc. But whatever the application, the final direction in which the observer sees the object can be expressed as direction cosines.

But where has this got us? Instead of two numbers – a longitude and a latitude – we now have three numbers to look after – the x-, y- and z-components – whose quadratic sum we have somehow to contrive to be unity. And, in addition to this apparent redundancy, most people find it harder to visualize problems in terms of $\left[x,y,z\right]$ than in $\left[\theta ,\varphi \right]$. Despite these objections, the vector approach turns out to have significant advantages over the spherical trigonometry approach:

• Vector formulae tend to be much more succinct; one vector operation is the equivalent of strings of sines and cosines.
• The formulae are as a rule rigorous, even at the poles.
• Accuracy is maintained all over the celestial sphere. When one Cartesian component is nearly unity and therefore insensitive to direction, the others become small and therefore more precise.
• Formulations usually deliver the quadrant of the result without the need for any inspection (except within the library function ATAN2).

A number of important transformations in positional astronomy turn out to be nothing more than changes of coordinate system, something which is especially convenient if the vector approach is used. A direction with respect to one triad can be expressed relative to another triad simply by multiplying the $\left[x,y,z\right]$ column vector by the appropriate $3\times 3$ orthogonal matrix (a tensor of Rank 2, or dyadic). The three rows of this rotation matrix are the vectors in the old coordinate system of the three new axes, and the transformation amounts to obtaining the dot-product of the direction-vector with each of the three new axes. Precession, nutation, $\left[h,\delta \right]$ to $\left[Az,El\right]$, $\left[\alpha ,\delta \right]$ to $\left[lII,bII\right]$ and so on are typical examples of the technique. A useful property of the rotation matrices is that they can be inverted simply by taking the transpose.

The elements of these vectors and matrices are assorted combinations of the sines and cosines of the various angles involved (hour angle, declination and so on, depending on which transformation is being applied). If you write out the matrix multiplications in full you get expressions which are essentially the same as the equivalent spherical trigonometry formulae. Indeed, many of the standard formulae of spherical trigonometry are most easily derived by expressing the problem initially in terms of vectors.

4.2.1 Using vectors

SLALIB provides conversions between spherical and vector form (sla_CS2C, sla_CC2S etc.), plus an assortment of standard vector and matrix operations (sla_VDV, sla_MXV etc.). There are also routines (sla_EULER etc.) for creating a rotation matrix from three Euler angles (successive rotations about specified Cartesian axes). Instead of Euler angles, a rotation matrix can be expressed as an axial vector (the pole of the rotation, and the amount of rotation), and routines are provided for this (sla_AV2M, sla_M2AV etc.).

Here is an example where spherical coordinates P1 and Q1 undergo a coordinate transformation and become P2 and Q2; the transformation consists of a rotation of the coordinate system through angles A, B and C about the z, new y and new z axes respectively:

Small adjustments to the direction of a position vector are often most conveniently described in terms of $\left[\Delta x,\Delta y,\Delta z\right]$. Adding the correction vector needs careful handling if the position vector is to remain of length unity, an advisable precaution which ensures that the $\left[x,y,z\right]$ components are always available to mean the cosines of the angles between the vector and the axis concerned. Two types of shifts are commonly used, the first where a small vector of arbitrary direction is added to the unit vector, and the second where there is a displacement in the latitude coordinate (declination, elevation etc.) alone.

For a shift produced by adding a small $\left[x,y,z\right]$ vector $$d to a unit vector ${}_{}$v1, the resulting vector ${}_{}$v2 has direction $<{}_{}$v1 + d > but is no longer of unit length. A better approximation is available if the result is multiplied by a scaling factor of $(1-$d ⋅v1), where the dot means scalar product. In Fortran:

The correction for diurnal aberration (discussed later) is an example of this form of shift.

As an example of the second kind of displacement we will apply a small change in elevation $\delta E$ to an $\left[Az,El\right]$ direction vector. The direction of the result can be obtained by making the allowable approximation $tan\delta E\approx \delta E$ and adding a adjustment vector of length $\delta E$ normal to the direction vector in the vertical plane containing the direction vector. The $z$-component of the adjustment vector is $\delta EcosE$, and the horizontal component is $\delta EsinE$ which has then to be resolved into $x$ and $y$ in proportion to their current sizes. To approximate a unit vector more closely, a correction factor of $cos\delta E$ can then be applied, which is nearly $\left(1-\delta E^2/2\right)$ for small $\delta E$. Expressed in Fortran, for initial vector V1X,V1Y,V1Z, change in elevation DEL (+ve $\equiv$ upwards), and result vector V2X,V2Y,V2Z:

An example of this type of shift is the correction for atmospheric refraction (see later). Depending on the relationship between $\delta E$ and $E$, special handling at the pole (the zenith for our example) may be required.

SLALIB includes routines for the case where both a position and a velocity are involved. The routines sla_CS2C6 and sla_CC62S convert from $\left[\theta ,\varphi ,\dot{\theta },\dot{\varphi }\right]$ to $\left[x,y,z,ẋ,ẏ,ż\right]$ and back; sla_DS2C6 and sla_DC62S are double precision equivalents.

4.3 Celestial Coordinate Systems

SLALIB has routines to perform transformations of celestial positions between different spherical coordinate systems, including those shown in the following table:

The routines sla_EQECL and sla_ECLEQ transform between ecliptic coordinates and $\left[\alpha ,\delta \right]$; there is also a routine for generating the equatorial to ecliptic rotation matrix for a given date: sla_ECMAT.

For conversion between Galactic coordinates and $\left[\alpha ,\delta \right]$ there are two sets of routines, depending on whether the $\left[\alpha ,\delta \right]$ is old-style, B1950, or new-style, J2000; sla_EG50 and sla_GE50 are $\left[\alpha ,\delta \right]$ to $\left[lII,bII\right]$ and vice versa for the B1950 case, while sla_EQGAL and sla_GALEQ are the J2000 equivalents.

Finally, the routines sla_GALSUP and sla_SUPGAL transform $\left[lII,bII\right]$ to de Vaucouleurs supergalactic longitude and latitude and vice versa.

It should be appreciated that the table, above, constitutes a gross oversimplification. Apparently simple concepts such as equator, equinox etc. are apt to be very hard to pin down precisely (polar motion, orbital perturbations …) and some have several interpretations, all subtly different. The various frames move in complicated ways with respect to one another or to the stars (themselves in motion). And in some instances the coordinate system is slightly distorted, so that the ordinary rules of spherical trigonometry no longer strictly apply.

These caveats apply particularly to the bewildering variety of different $\left[\alpha ,\delta \right]$ systems that are in use. Figure 1 shows how some of these systems are related, to one another and to the direction in which a celestial source actually appears in the sky. At the top of the diagram are the various sorts of mean place found in star catalogues and papers;² at the bottom is the observed $\left[Az,El\right]$, where a perfect theodolite would be pointed to see the source; and in the body of the diagram are the intermediate processing steps and coordinate systems. To help understand this diagram, and the SLALIB routines that can be used to carry out the various calculations, we will look at the coordinate systems involved, and the astronomical phenomena that affect them.

4.4 Precession and Nutation

Right ascension and declination, ($\left[\alpha ,\delta \right]$), are the names of the longitude and latitude in a spherical polar coordinate system based on the Earth’s axis of rotation. The zero point of $\alpha$ is the point of intersection of the celestial equator and the ecliptic (the apparent path of the Sun through the year) where the Sun moves into the northern hemisphere. This point is called the first point of Aries, the vernal equinox (with apologies to southern-hemisphere readers) or simply the equinox.³

This simple picture is unfortunately complicated by the difficulty of defining a suitable equator and equinox. One problem is that the Sun’s apparent diurnal and annual motions are not completely regular, due to the ellipticity of the Earth’s orbit and its continuous disturbance by the Moon and planets. This is dealt with by separating the motion into (i) a smooth and steady mean Sun and (ii) a set of periodic corrections and perturbations; only the former is involved in establishing reference systems and time scales. A second, far larger problem, is that the celestial equator and the ecliptic are both moving with respect to the stars. These motions arise because of the gravitational interactions between the Earth and the other solar-system bodies.

By far the largest effect is the so-called “precession of the equinoxes”, where the Earth’s rotation axis sweeps out a cone centred on the ecliptic pole, completing one revolution in about 26,000 years. The cause of the motion is the torque exerted on the distorted and spinning Earth by the Sun and the Moon. Consider the effect of the Sun alone, at or near the northern summer solstice. The Sun ‘sees’ the top (north pole) of the Earth tilted towards it (by about $23\circ .5$, the obliquity of the ecliptic), and sees the nearer part of the Earth’s equatorial bulge below centre and the further part above centre. Although the Earth is in free fall, the gravitational force on the nearer part of the equatorial bulge is greater than that on the further part, and so there is a net torque acting as if to eliminate the tilt. Six months later the same thing is happening in reverse, except that the torque is still trying to eliminate the tilt. In between (at the equinoxes) the torque shrinks to zero. A torque acting on a spinning body is gyroscopically translated into a precessional motion of the spin axis at right-angles to the torque, and this happens to the Earth. The motion varies during the year, going through two maxima, but always acts in the same direction. The Moon produces the same effect, adding a contribution to the precession which peaks twice per month. The Moon’s proximity to the Earth more than compensates for its smaller mass and gravitational attraction, so that it in fact contributes most of the precessional effect.

The complex interactions between the three bodies produce a precessional motion that is wobbly rather than completely smooth. However, the main 26,000-year component is on such a grand scale that it dwarfs the remaining terms, the biggest of which has an amplitude of only $11$${}^{\prime }{}^{\prime }$ and a period of about 18.6 years. This difference of scale makes it convenient to treat these two components of the motion separately. The main 26,000-year effect is called luni-solar precession; the smaller, faster, periodic terms are called the nutation.

Note that precession and nutation are simply different frequency components of the same physical effect. It is a common misconception that precession is caused by the Sun and nutation is caused by the Moon. In fact the Moon is responsible for two-thirds of the precession, and, while it is true that much of the complex detail of the nutation is a reflection of the intricacies of the lunar orbit, there are nonetheless important solar terms in the nutation.

In addition to and quite separate from the precession-nutation effect, the orbit of the Earth-Moon system is not fixed in orientation, a result of the attractions of the planets. This slow (about ′′05 per year) secular rotation of the ecliptic about a slowly-moving diameter is called, confusingly, planetary precession and, along with the luni-solar precession is included in the general precession. The equator and ecliptic as affected by general precession are what define the various “mean” $\left[\alpha ,\delta \right]$ reference frames.

The models for precession and nutation come from a combination of observation and theory, and are subject to continuous refinement. Nutation models in particular have reached a high degree of sophistication, taking into account such things as the non-rigidity of the Earth and the effects of the planets; SLALIB’s nutation model (SF2001) involves 194 terms in each of $\psi$ (longitude) and $ϵ$ (obliquity), some as small as a few microarcseconds.

4.4.1 SLALIB support for precession and nutation

SLALIB offers a choice of three precession models:

• The old Bessel-Newcomb, pre IAU 1976, “FK4” model, used for B1950 star positions and other pre-1984.0 purposes (sla_PREBN).
• The new Fricke, IAU 1976, “FK5” model, used for J2000 star positions and other post-1984.0 purposes (sla_PREC).
• A model published by Simon et al. which is more accurate than the IAU 1976 model and which is suitable for long periods of time (sla_PRECL).

In each case, the named SLALIB routine generates the $\left(3\times 3\right)$ precession matrix for a given start and finish time. For example, here is the Fortran code for generating the rotation matrix which describes the precession between the epochs J2000 and J1985.372 (IAU 1976 model):

It is instructive to examine the resulting matrix:

Note that the diagonal elements are close to unity, and the other elements are small. This shows that over an interval as short as 15 years the precession isn’t going to move a position vector very far (in this case about $0\circ .2$).

For convenience, a direct $\left[\alpha ,\delta \right]$ to $\left[\alpha ,\delta \right]$ precession routine is also provided (sla_PRECES), suitable for either the old or the new system (but not a mixture of the two).

SLALIB provides two nutation models, the old IAU 1980 model, implemented in the routine slaNutc80, and a much more accurate newer theory, SF2001, implemented in the routine slaNutc. Both return the components of nutation in longitude and latitude (and also provide the obliquity) from which a nutation matrix can be generated by calling slaDeuler (and from which the equation of the equinoxes, described later, can be found). Alternatively, the SF2001 nutation matrix can be generated in a single call by using slaNut. The SF2001 nutation theory includes components that correct for errors in the IAU 1976 precession and also for the $\sim 23$mas displacement between the mean J2000 and ICRS coordinate systems, achieving a final accuracy well under 1 mas in the present era.

A rotation matrix for applying the entire precession-nutation transformation in one go can be generated by calling sla_PRENUT.

4.5 Mean Places

From a classical standpoint, the main effect of the precession-nutation is an increase of about $50$${}^{\prime }{}^{\prime }$/year in the ecliptic longitudes of the stars. It is therefore essential, when reporting the position of an astronomical target, to qualify the coordinates with a date, or epoch. Specifying the epoch ties down the equator and equinox which define the $\left[\alpha ,\delta \right]$ coordinate system that is being used. ⁴ For simplicity, only the smooth and steady “precession” part of the complete precession-nutation effect is included, thereby defining what is called the mean equator and equinox for the epoch concerned. We say a star has a mean place of (for example) $12h07m$$58$$s$$.09$ $-19\circ 44^{\prime }$′′371 “with respect to the mean equator and equinox of epoch J2000”. The short way of saying this is “$\left[\alpha ,\delta \right]$ equinox J2000” (not “$\left[\alpha ,\delta \right]$ epoch J2000”, which means something different to do with proper motion).

4.6 Epoch

The word “epoch” just means a significant moment in time, and can be supplied in a variety of forms, using different calendar systems and time scales.

For the purpose of specifying the epochs associated with the mean place of a star, two conventions exist. Both sorts of epoch superficially resemble years AD but are not tied to the civil (Gregorian) calendar; to distinguish them from ordinary calendar-years there is often a “.0” suffix (as in “1950.0”), although any other fractional part is perfectly legal (e.g. 1987.5).

The older system, Besselian epoch, is defined in such a way that its units are tropical years of about 365.2422 days and its time scale is the obsolete Ephemeris Time. The start of the Besselian year is the moment when the ecliptic longitude of the mean Sun is $280^{\circ }$; this happens near the start of the calendar year (which is why $280^{\circ }$ was chosen).

The new system, Julian epoch, was adopted as part of the IAU 1976 revisions (about which more will be said in due course) and came formally into use at the beginning of 1984. It uses the Julian year of exactly 365.25 days; Julian epoch 2000 is defined to be 2000 January 1.5 in the TT time scale.

For specifying mean places, various standard epochs are in use, the most common ones being Besselian epoch 1950.0 and Julian epoch 2000.0. To distinguish the two systems, Besselian epochs are now prefixed “B” and Julian epochs are prefixed “J”. Epochs without an initial letter can be assumed to be Besselian if before 1984.0, otherwise Julian. These details are supported by the SLALIB routines sla_DBJIN (decodes numbers from a character string, accepting an optional leading B or J), sla_KBJ (decides whether B or J depending on prefix or range) and sla_EPCO (converts one epoch to match another).

SLALIB has four routines for converting Besselian and Julian epochs into other forms. The functions sla_EPB2D and sla_EPJ2D convert Besselian and Julian epochs into MJD; the functions sla_EPB and sla_EPJ do the reverse. For example, to express B1950 as a Julian epoch:

(The answer is J1949.99979.)

4.7 Proper Motion

Stars in catalogues usually have, in addition to the $\left[\alpha ,\delta \right]$ coordinates, a proper motion $\left[{\mu }_{\alpha },{\mu }_{\delta }\right]$. This is an intrinsic motion of the star across the background. Very few stars have a proper motion which exceeds $1$${}^{\prime }{}^{\prime }$/year, and most are far below this level. A star observed as part of normal astronomy research will, as a rule, have a proper motion which is unknown.

Mean $\left[\alpha ,\delta \right]$ and rate of change are not sufficient to pin down a star; the epoch at which the $\left[\alpha ,\delta \right]$ was or will be correct is also needed. Note the distinction between the epoch which specifies the coordinate system and the epoch at which the star passed through the given $\left[\alpha ,\delta \right]$. The full specification for a star is $\left[\alpha ,\delta \right]$, proper motions, equinox and epoch (plus something to identify which set of models for the precession etc. is being used – see the next section). For convenience, coordinates given in star catalogues are almost always adjusted to make the equinox and epoch the same – for example B1950 in the case of the SAO catalogue.

SLALIB provides one routine to handle proper motion on its own, sla_PM. Proper motion is also allowed for in various other routines as appropriate, for example sla_MAP and sla_FK425. Note that in all SLALIB routines which involve proper motion the units are radians per year and the $\alpha$ component is in the form $\dot{\alpha }$ (i.e. big numbers near the poles). Some star catalogues have proper motion per century, and in some catalogues the $\alpha$ component is in the form $\dot{\alpha }cos\delta$ (i.e. angle on the sky).

4.8 Parallax and Radial Velocity

For the utmost accuracy and the nearest stars, allowance can be made for annual parallax and for the effects of perspective on the proper motion.

Parallax is appreciable only for nearby stars; even the nearest, Proxima Centauri, is displaced from its average position by less than an arcsecond as the Earth revolves in its orbit.

For stars with a known parallax, knowledge of the radial velocity allows the proper motion to be expressed as an actual space motion in 3 dimensions. The proper motion is, in fact, a snapshot of the transverse component of the space motion, and in the case of nearby stars will change with time due to perspective.

SLALIB does not provide facilities for handling parallax and radial-velocity on their own, but their contribution is allowed for in such routines as sla_PM, sla_MAP and sla_FK425. Catalogue mean places do not include the effects of parallax and are therefore barycentric; when pointing telescopes etc. it is usually most efficient to apply the slowly-changing parallax correction to the mean place of the target early on and to work with the geocentric mean place. This latter approach is implied in Figure 1.

4.9 Aberration

The finite speed of light combined with the motion of the observer around the Sun during the year causes apparent displacements of the positions of the stars. The effect is called the annual aberration (or “stellar” aberration). Its maximum size, about ′′205, occurs for stars $90^{\circ }$ from the point towards which the Earth is headed as it orbits the Sun; a star exactly in line with the Earth’s motion is not displaced. To receive the light of a star, the telescope has to be offset slightly in the direction of the Earth’s motion. A familiar analogy is the need to tilt your umbrella forward when on the move, to avoid getting wet. This classical model is, in fact, misleading in the context of light as opposed to rain, but happens to give the same answer as a relativistic treatment to first order (better than 1 milliarcsecond).

Before the IAU 1976 resolutions, different values for the approximately ′′205 aberration constant were employed at different times, and this can complicate comparisons between different catalogues. Another complication comes from the so-called E-terms of aberration, that small part of the annual aberration correction that is a function of the eccentricity of the Earth’s orbit. The E-terms, maximum amplitude about ′′03, happen to be approximately constant for a given star, and so they used to be incorporated in the catalogue $\left[\alpha ,\delta \right]$ to reduce the labour of converting to and from apparent place. The E-terms can be removed from a catalogue $\left[\alpha ,\delta \right]$ by calling sla_SUBET or applied (for example to allow a pulsar timing-position to be plotted on a B1950 finding chart) by calling sla_ADDET; the E-terms vector itself can be obtained by calling sla_ETRMS. Star positions post IAU 1976 are free of these distortions, and to apply corrections for annual aberration involves the actual barycentric velocity of the Earth rather than the use of canonical circular-orbit models.

The annual aberration is the aberration correction for an imaginary observer at the Earth’s centre. The motion of a real observer around the Earth’s rotation axis in the course of the day makes a small extra contribution to the total aberration effect called the diurnal aberration. Its maximum amplitude is about ′′03.

No SLALIB routine is provided for calculating the aberration on its own, though the required velocity vectors can be generated using sla_EVP (or sla_EPV) and sla_GEOC. Annual and diurnal aberration are allowed for where required, for example in sla_MAP etc. and sla_AOP etc. Note that this sort of aberration is different from the planetary aberration, which is the apparent displacement of a solar-system body, with respect to the ephemeris position, as a consequence of the motion of both the Earth and the source. The planetary aberration can be computed either by correcting the position of the solar-system body for light-time, followed by the ordinary stellar aberration correction, or more directly by expressing the position and velocity of the source in the observer’s frame and correcting for light-time alone.

4.10 Different Sorts of Mean Place

A confusing aspect of the mean places used in the pre-ICRS era is that they are sensitive to the precise way they were determined. A mean place is not directly observable, even with fundamental instruments such as transit circles, and to produce one will involve relying on some existing star catalogue, for example the fundamental catalogues FK4 and FK5, and applying given mathematical models of precession, nutation, aberration and so on. Note in particular that no star catalogue, even a fundamental catalogue such as FK4 or FK5, defines a coordinate system, strictly speaking; it is merely a list of star positions and proper motions. However, once the stars from a given catalogue are used as position calibrators, e.g. for transit-circle observations or for plate reductions, then a broader sense of there being a coordinate grid naturally arises, and such phrases as “in the system of the FK4” can legitimately be employed. However, there is no formal link between the two concepts – no “standard least squares fit” between reality and the inevitably flawed catalogues. All such catalogues suffer at some level from systematic, zonal distortions of both the star positions and of the proper motions, and include measurement errors peculiar to individual stars.

Many of these complications are of little significance except to specialists. However, observational astronomers cannot escape exposure to at least the two main varieties of mean place, loosely called FK4 and FK5, and should be aware of certain pitfalls. For most practical purposes the more recent system, FK5, is free of surprises and tolerates naive use well. FK4, in contrast, contains two important traps:

• The FK4 system rotates at about ′′05 per century relative to distant galaxies. This is manifested as a systematic distortion in the proper motions of all FK4-derived catalogues, which will in turn pollute any astrometry done using those catalogues. For example, FK4-based astrometry of a QSO using plates taken decades apart will reveal a non-zero fictitious proper motion, and any FK4 star which happens to have zero proper motion is, in fact, slowly moving against the distant background. The FK4 frame rotates because it was established before the nature of the Milky Way, and hence the existence of systematic motions of nearby stars, had been recognized.
• Star positions in the FK4 system are part-corrected for annual aberration (see above) and embody the so-called E-terms of aberration.

The change from the old FK4-based system to FK5 occurred at the beginning of 1984 as part of a package of resolutions made by the IAU in 1976, along with the adoption of J2000 as the reference epoch. Star positions in the newer, FK5, system are free from the E-terms, and the system is a much better approximation to an inertial frame – about five times better (and ICRS is hundreds of times better still).

It may occasionally be convenient to specify the FK4 fictitious proper motion directly. In FK4, the centennial proper motion of (for example) a QSO is:

${\mu }_{\alpha }=-$$0$$s$$.015869$$+(($$0$$s$$.029032$$sin\alpha +$$0$$s$$.000340$$cos\alpha )sin\delta -$$0$$s$$.000105$$cos\alpha -$$0$$s$$.000083$$sin\alpha )sec\delta$
${\mu }_{\delta }=+$′′043549$cos\alpha -$′′000510$sin\alpha +($′′000158$sin\alpha -$′′000125$cos\alpha )sin\delta -$′′000066$cos\delta$

4.11 Mean Place Transformations

Figure 1 is based upon three varieties of mean $\left[\alpha ,\delta \right]$ all of which are of practical significance to observing astronomers in the present era:

• Old style (FK4) with known proper motion in the FK4 system, and with parallax and radial velocity either known or assumed zero.
• Old style (FK4) with zero proper motion in FK5, and with parallax and radial velocity assumed zero.
• New style (FK5 or, loosely, ICRS) with proper motion, parallax and radial velocity either known or assumed zero.

The figure outlines the steps required to convert positions in any of these systems to a J2000 $\left[\alpha ,\delta \right]$ for the current epoch, as might be required in a telescope-control program for example. Most of the steps can be carried out by calling a single SLALIB routine; there are other SLALIB routines which offer set-piece end-to-end transformation routines for common cases. Note, however, that SLALIB does not set out to provide the capability for arbitrary transformations of star-catalogue data between all possible systems of mean $\left[\alpha ,\delta \right]$. Only in the (common) cases of FK4, equinox and epoch B1950, to FK5, equinox and epoch J2000, and vice versa are proper motion, parallax and radial velocity transformed along with the star position itself, the focus of SLALIB support.

As an example of using SLALIB to transform mean places, here is Fortran code that implements the top-left path of Figure 1. An FK4 $\left[\alpha ,\delta \right]$ of arbitrary equinox and epoch and with known proper motion and parallax is transformed into an FK5 J2000 $\left[\alpha ,\delta \right]$ for the current epoch. As a test star we will use $\alpha =$$16h09m$$55$$s$$.13$, $\delta =$$-75\circ 59^{\prime }$′′272, equinox 1900, epoch 1963.087, ${\mu }_{\alpha }=$$-0$$s$$.0312$$/y$, ${\mu }_{\delta }=$′′+0103$/y$, parallax = ′′0062, radial velocity = $-34.22$ km/s. The date of observation is 1994.35.

It is interesting to look at how the $\left[\alpha ,\delta \right]$ changes during the course of the calculation:

Other remarks about the above (unusually complicated) example:

• If the original equinox and epoch were B1950, as is quite likely, then it would be unnecessary to treat space motions and E-terms explicitly. Transformation to FK5 J2000 could be accomplished simply by calling sla_FK425, after which a call to sla_PM and the parallax code would complete the work.
• The rigorous treatment of the E-terms has only a small effect on the result. Such refinements are, nevertheless, worthwhile in order to facilitate comparisons and to increase the chances that star positions from different suppliers are compatible.
• The FK4 to FK5 transformations, sla_FK425 and sla_FK45Z, are not as is sometimes assumed simply 50 years of precession, though this indeed accounts for most of the change. The transformations also include adjustments to the equinox, a revised precession model, elimination of the E-terms, a change to the proper-motion time unit and so on. The reason there are two routines rather than just one is that the FK4 frame rotates relative to the background, whereas the FK5 frame is a much better approximation to an inertial frame, and zero proper motion in FK4 does not, therefore, mean zero proper motion in FK5. SLALIB also provides two routines, sla_FK524 and sla_FK54Z, to perform the inverse transformations.
• Some star catalogues (FK4 itself is one) were constructed using slightly different procedures for the polar regions compared with elsewhere. SLALIB ignores this inhomogeneity and always applies the standard transformations, irrespective of location on the celestial sphere.

4.12 Mean Place to Apparent Place

The geocentric apparent place of a source, or apparent place for short, is the $\left[\alpha ,\delta \right]$ if viewed from the centre of the Earth, with respect to the true equator and equinox of date. Transformation of an FK5 mean $\left[\alpha ,\delta \right]$, equinox J2000, current epoch, to apparent place involves the following effects:

• Light deflection – the gravitational lens effect of the sun.
• Annual aberration.
• Precession-nutation.

The light deflection is seldom significant. Its value at the limb of the Sun is about ′′174; it falls off rapidly with distance from the Sun and has shrunk to about ′′002 at an elongation of $20^{\circ }$.

As already described, the annual aberration is a function of the Earth’s velocity relative to the solar system barycentre (available through the SLALIB routines sla_EVP and sla_EPV) and produces shifts of up to about ′′205.

The precession-nutation, from J2000 to the current epoch, is expressed by a rotation matrix which is available through the SLALIB routine sla_PRENUT.

The whole mean-to-apparent transformation can be done using the SLALIB routine sla_MAP. As a demonstration, here is a program which lists the North Polar Distance ($90^{\circ }-\delta$) of Polaris for the decade of closest approach to the Pole:

For cases where the transformation has to be repeated for different times or for more than one star, the straightforward sla_MAP approach is apt to be wasteful as both the Earth velocity and the precession-nutation matrix can be re-calculated relatively infrequently without ill effect. A more efficient method is to perform the target-independent calculations only when necessary, by calling sla_MAPPA, and then to use either sla_MAPQKZ, when only the $\left[\alpha ,\delta \right]$ is known, or sla_MAPQK, when full catalogue positions, including proper motion, parallax and radial velocity, are available. How frequently to call sla_MAPPA depends on the accuracy objectives; once per night will deliver sub-arcsecond accuracy for example.

The routines sla_AMP and sla_AMPQK allow the reverse transformation, from apparent to mean place.

4.13 Apparent Place to Observed Place

The observed place of a source is its position as seen by a perfect theodolite at the location of the observer. Transformation of an apparent $\left[\alpha ,\delta \right]$ to observed place involves the following effects:

• $\left[\alpha ,\delta \right]$ to $\left[h,\delta \right]$.
• Diurnal aberration.
• $\left[h,\delta \right]$ to $\left[Az,El\right]$.
• Refraction.

The transformation from apparent $\left[\alpha ,\delta \right]$ to apparent $\left[h,\delta \right]$ is made by allowing for Earth rotation through the sidereal time, $\theta$:

For this equation to work, $\alpha$ must be the apparent right ascension for the time of observation, and $\theta$ must be the local apparent sidereal time. The latter is obtained as follows:

(1)

from civil time obtain the coordinated universal time, UTC (more later on this);

(2)

add the UT1$-$UTC (typically a few tenths of a second) to give the UT;

(3)

from the UT compute the Greenwich mean sidereal time (using sla_GMST);

(4)

add the observer’s (east) longitude, giving the local mean sidereal time;

(5)

add the equation of the equinoxes (using sla_EQEQX).

The equation of the equinoxes ($=\Delta \psi cosϵ$ plus small terms) is the effect of nutation on the sidereal time. Its value is typically a second or less. It is interesting to note that if the object of the exercise is to transform a mean place all the way into an observed place (very often the case), then the equation of the equinoxes and the longitude component of nutation can both be omitted, removing a great deal of computation. However, SLALIB follows the normal convention and works via the apparent place.

Note that for very precise work the observer’s longitude should be corrected for polar motion. This can be done with sla_POLMO. The corrections are always less than about ′′03, and are futile unless the position of the observer’s telescope is known to better than a few metres.

Tables of observed and predicted UT1$-$UTC corrections and polar motion data are published every few weeks by the International Earth Rotation Service.

The transformation from apparent $\left[h,\delta \right]$ to topocentric $\left[h,\delta \right]$ consists of allowing for diurnal aberration. This effect, maximum amplitude ′′02, was described earlier. There is no specific SLALIB routine for computing the diurnal aberration, though the routines sla_AOP etc. include it, and the required velocity vector can be determined by calling sla_GEOC.

The next stage is the major coordinate rotation from local equatorial coordinates $\left[h,\delta \right]$ into horizon coordinates. The SLALIB routines sla_E2H etc. can be used for this. For high-precision applications the mean geodetic latitude should be corrected for polar motion.

4.13.1 Refraction

The final correction is for atmospheric refraction. This effect, which depends on local meteorological conditions and the effective colour of the source/detector combination, increases the observed elevation of the source by a significant effect even at moderate zenith distances, and near the horizon by over $0\circ .5$. The amount of refraction can by computed by calling the SLALIB routine sla_REFRO; however, this requires as input the observed zenith distance, which is what we are trying to predict. For high precision it is therefore necessary to iterate, using the topocentric zenith distance as the initial estimate of the observed zenith distance.

The full sla_REFRO refraction calculation is onerous, and for zenith distances of less than, say, $75^{\circ }$ the following model can be used instead:

where ${\zeta }_{vac}$ is the topocentric zenith distance (i.e. in vacuo), ${\zeta }_{obs}$ is the observed zenith distance (i.e. affected by refraction), and $A$ and $B$ are constants, about $60$${}^{\prime }{}^{\prime }$ and ′′-006 respectively for a sea-level site. The two constants can be calculated for a given set of conditions by calling either sla_REFCO or sla_REFCOQ.

sla_REFCO works by calling sla_REFRO for two zenith distances and fitting $A$ and $B$ to match. The calculation is onerous, but delivers accurate results whatever the conditions. sla_REFCOQ uses a direct formulation of $A$ and $B$ and is much faster; it is slightly less accurate than sla_REFCO but more than adequate for most practical purposes.

Like the full refraction model, the two-term formulation works in the wrong direction for our purposes, predicting the in vacuo (topocentric) zenith distance given the refracted (observed) zenith distance, rather than vice versa. The obvious approach of interchanging ${\zeta }_{vac}$ and ${\zeta }_{obs}$ and reversing the signs, though approximately correct, gives avoidable errors which are just significant in some applications; for example about ′′02 at $70^{\circ }$ zenith distance. A much better result can easily be obtained, by using one Newton-Raphson iteration as follows:

The effect of refraction can be applied to an unrefracted zenith distance by calling sla_REFZ or to an unrefracted $\left[x,y,z\right]$ by calling sla_REFV. Over most of the sky these two routines deliver almost identical results, but beyond $\zeta =83^{\circ }$ sla_REFV becomes unacceptably inaccurate while sla_REFZ remains usable. (However sla_REFV is significantly faster, which may be important in some applications.) SLALIB also provides a routine for computing the airmass, the function sla_AIRMAS.

The refraction “constants” returned by sla_REFCO and sla_REFCOQ are slightly affected by colour, especially at the blue end of the spectrum. Where values for more than one wavelength are needed, rather than calling sla_REFCO several times it is more efficient to call sla_REFCO just once, for a selected “base” wavelength, and then to call sla_ATMDSP once for each wavelength of interest.

All the SLALIB refraction routines work for radio wavelengths as well as the optical/IR band. The radio refraction is very dependent on humidity, and an accurate value must be supplied. There is no wavelength dependence, however. The choice of optical/IR or radio is made by specifying a wavelength greater than $100\mu m$ for the radio case.

4.13.2 Efficiency considerations

The complete apparent place to observed place transformation can be carried out by calling sla_AOP. For improved efficiency in cases of more than one star or a sequence of times, the target-independent calculations can be done once by calling sla_AOPPA, the time can be updated by calling sla_AOPPAT, and sla_AOPQK can then be used to perform the apparent-to-observed transformation. The reverse transformation is available through sla_OAP and sla_OAPQK. (n.b. These routines use accurate but computationally-expensive refraction algorithms for zenith distances beyond about $76^{\circ }$. For many purposes, in-line code tailored to the accuracy requirements of the application will be preferable, for example ignoring polar motion, omitting diurnal aberration and using sla_REFZ to apply the refraction.)

4.14 The Hipparcos Catalogue and the ICRS

With effect from the beginning of 1998, the IAU adopted a new reference system to replace FK5 J2000. The new system, called the International Celestial Reference System (ICRS), differs profoundly from all predecessors in that the link with solar-system dynamics was broken; the ICRS axes are defined in terms of the coordinates of a set of extragalactic sources, not in terms of the mean equator and equinox at a given reference epoch. Although the ICRS and FK5 coordinates of any given object are almost the same, the orientation of the new frame was essentially arbitrary, and the close match to FK5 J2000 was contrived purely for reasons of continuity and convenience.

A distinction is made between the reference system (the ICRS) and frame (ICRF). The ICRS is the set of prescriptions and conventions together with the modelling required to define, at any time, a triad of axes. The ICRF is a practical realization, and currently consists of a catalogue of equatorial coordinates for 608 extragalactic radio sources observed by VLBI.

The best optical realization of the ICRF currently available is the Hipparcos catalogue. The extragalactic sources were not directly observable by the Hipparcos satellite and so the link from Hipparcos to ICRF was established through a variety of indirect techniques: VLBI and conventional interferometry of radio stars, photographic astrometry and so on. The Hipparcos frame is aligned to the ICRF to within about 0.5 mas and 0.5 mas/year (at epoch 1991.25).

The Hipparcos catalogue includes all of the FK5 stars, which has enabled the orientation and spin of the latter to be studied. At epoch J2000, the misalignment of the FK5 frame with respect to Hipparcos (and hence ICRS) are about 32 mas and 1 mas/year respectively. Consequently, for many practical purposes, including pointing telescopes, the IAU 1976-1982 conventions on reference frames and Earth orientation remain adequate and there is no need to change to Hipparcos coordinates, new precession-nutation models and so on. However, for the most exacting astrometric applications, SLALIB provides some support for Hipparcos coordinates in the form of four new routines: sla_FK52H and sla_H2FK5, which transform FK5 positions and proper motions to the Hipparcos frame and vice versa, and sla_FK5HZ and sla_HFK5Z, where the transformations are for stars whose Hipparcos proper motion is zero.

Further information on the ICRS can be found in the paper by M. Feissel and F. Mignard, Astron. Astrophys. 331, L33-L36 (1988).

4.15 Time Scales

SLALIB provides for transformation between several time scales, and involves use of one or two others. The full list is as follows:

• TAI: International Atomic Time
• UTC: Coordinated Universal Time
• TT: Terrestrial Time
• TDB: Barycentric Dynamical Time.
• UT: Universal Time
• GMST: Greenwich mean sidereal time
• GAST (or GST): Greenwich apparent sidereal time.
• LAST: local apparent sidereal time

Strictly speaking, UT and the sidereal times are not times in the physics sense, but angles that describe Earth rotation.

Three obsolete time scales should be mentioned here to avoid confusion.

• GMT: Greenwich Mean Time – can mean either UTC or UT.
• ET: Ephemeris Time – more or less the same as either TT or TDB.
• TDT: Terrestrial Dynamical Time – former name of TT.

time scales that have no SLALIB support at present:

• Any form of local civil time (BST, PDT etc.)
• TCG: geocentric coordinate time.
• TCB: barycentric coordinate time.

4.15.1 Atomic Time: TAI

International Atomic Time, TAI, is a “laboratory” time scale with no link to astronomical observations except in an historical sense. Its unit is the SI second, which is defined in terms of a specific number of wavelengths of the radiation produced by a certain electronic transition in the caesium 133 atom. It is realized through a changing population of high-precision atomic clocks held at standards institutes in various countries. There is an elaborate process of continuous intercomparison, leading to a weighted average of all the clocks involved.

Though TAI shares the same second as the more familiar UTC, the two time scales are noticeably separated in epoch because of the build-up of leap seconds (see the next section). At the time of writing, UTC lags over half a minute behind TAI.

For any given date, the difference TAI$-$UTC can be obtained by calling the SLALIB function sla_DAT. Note, however, that an up-to-date copy of the function must be used if the most recent leap seconds are required. For applications where this is critical, mechanisms independent of SLALIB and under local control must be set up; in such cases sla_DAT can be useful as an independent check, for test dates within the range of the available version. Up-to-date information on TAI$-$UTC is available from ftp://maia.usno.navy.mil/ser7/tai-utc.dat.

4.15.2 Universal Time: UT, UTC

Universal Time, UT, or more specifically UT1, is in effect mean solar time and is really an expression of Earth rotation rather than a measure of time. Originally defined in terms of a point in the sky called “the fictitious mean Sun”, UT is now defined through its relationship with Earth rotation angle (formerly sidereal time). Because the Earth’s rotation rate is slightly irregular and gradually decreasing,⁵ the UT second is not precisely matched to the SI second. This makes UT itself unsuitable for use as a time scale.

That role is instead taken by Coordinated Universal Time, UTC, which is clock-based and is the foundation of civil timekeeping. Most time zones differ from UTC by an integer number of hours, though a few (e.g. parts of Canada and Australia) differ by $n+0.5$ hours. Since its introduction, UTC has been kept roughly in step with UT by a variety of adjustments that are agreed in advance and then carried out in a coordinated manner by the timekeeping communities of different countries—hence the name. Though rate changes were used in the past, nowadays all such adjustments are made by occasionally inserting a whole second. This procedure is called a leap second. Because the day length is now slightly longer than 86400 SI seconds, a leap second amounts to stopping the UTC clock for a second to let the Earth catch up.

You need UT1 in order to point a telescope or antenna at a celestial target. To obtain it starting from UTC, you have to look up the value of UT1$-$UTC for the date concerned in tables published by the International Earth Rotation and reference frames Service; this quantity, kept in the range $\pm$$0$$s$$.9$ by means of leap seconds, is then added to the UTC. The quantity UT1$-$UTC, which typically changes by of order 1 ms per day, can be obtained only by observation (VLBI using extragalactic radio sources), though seasonal trends are well known and the IERS listings are able to predict some way into the future with adequate accuracy for pointing telescopes.

UTC leap seconds are introduced as necessary, usually at the end of December or June. Because on the average the solar day is slightly longer than the nominal 86,400 SI seconds, leap seconds are always positive; however, provision exists for negative leap seconds if needed. The form of a leap second can be seen from the following description of the end of June 1994:

Note that UTC has to be expressed as hours, minutes and seconds (or at least in seconds for a given date) if leap seconds are to be taken into account in the correct manner. It is improper to express a UTC as a Julian Date, for example, because there will be an ambiguity during a leap second (in the above example, 1994 June 30 $23h59m$$60$$s$$.0$ and 1994 July 1 $00h00m$$00$$s$$.0$ would both come out as MJD 49534.00000). Although in the vast majority of cases this won’t matter, there are potential problems in on-line data acquisition systems and in applications involving taking the difference between two times. Note that although the functions sla_DAT and sla_DTT expect UTC in the form of an MJD, the meaning here is really a whole-number date rather than a time. Though the functions will accept a fractional part and will almost always function correctly, on a day which ends with a leap second incorrect results would be obtained during the leap second itself because by then the MJD would have moved into the next day.

4.15.3 Sidereal Time: GMST, LAST etc.

Sidereal time is like the time of day but relative to the stars rather than to the Sun. After one sidereal day the stars come back to the same place in the sky, apart from sub-arcsecond precession effects. Because the Earth rotates faster relative to the stars than to the Sun by one day per year, the sidereal second is shorter than the solar second; the ratio is about 0.9973.

The Greenwich mean sidereal time, GMST, is linked to UT1 by a numerical formula which is implemented in the SLALIB functions sla_GMST and sla_GMSTA. There are, of course, no leap seconds in GMST, but the sidereal second (measured in SI seconds) changes in length along with the UT1 second, and also varies over long periods of time because of slow changes in the Earth’s orbit. This makes sidereal time unsuitable for everything except predicting the apparent directions of celestial sources, in other words as an angle rather than a time.

The local apparent sidereal time, LAST, is the apparent right ascension of the local meridian, from which the hour angle of any star can be determined knowing its right ascension. LAST can be obtained from the GMST by adding the east longitude (corrected for polar motion in precise work) and the equation of the equinoxes. The latter, already described, is an aspect of the nutation effect and can be predicted by calling the SLALIB function sla_EQEQX or, neglecting certain very small terms, by calling sla_NUTC and using the expression $\Delta \psi cosϵ$.

GAST, or plain GST, is GMST plus the equation of the equinoxes.

4.15.4 Dynamical Time: TT, TDB

Dynamical time (formerly Ephemeris Time, ET) is the independent variable in the theories which describe the motions of bodies in the solar system. When using published formulae or tables that model the position of the Earth in its orbit, for example, or look up the Moon’s position in a precomputed ephemeris, the date and time must be in terms of one of the dynamical time scales. It is a common but understandable mistake to use UTC directly, in which case the results will be over a minute out (at the time of writing).

It is not hard to see why such time scales are necessary. UTC would clearly be unsuitable as the argument of an ephemeris because of leap seconds. A solar-system ephemeris based on UT1 or sidereal time would somehow have to include the unpredictable variations of the Earth’s rotation. TAI would work, but in principle the ephemeris and the ensemble of atomic clocks would eventually drift apart. In effect, the ephemeris is a clock, with the bodies of the solar system the hands from which the ephemeris time is read.

Only two of the dynamical time scales are of any great importance to observational astronomers, TT and TDB.

Terrestrial Time, TT, is the theoretical time scale of apparent geocentric ephemerides of solar system bodies. It applies to clocks at sea-level, and for practical purposes it is tied to Atomic Time TAI through the formula TT $=$ TAI $+$ $32$$s$$.184$. In practice, therefore, the units of TT are ordinary SI seconds, and the offset of $32$$s$$.184$ with respect to TAI is fixed. The SLALIB function sla_DTT returns TT$-$UTC for a given UTC (n.b. sla_DTT calls sla_DTT, and the latter must be an up-to-date version if recent leap seconds are to be taken into account).

Barycentric Dynamical Time, TDB, is a coordinate time, suitable for labelling events that are most simply described in a context where the bodies of the solar system are absent. Applications include the emission of pulsar radiation and the motions of the solar-system bodies themselves. When the readings of the observer’s TT clock are labelled using such a coordinate time, differences are seen because the clock is affected by its speed in the barycentric coordinate system and the gravitational potential in which it is immersed. Equivalently, observations of pulsars expressed in TT would display similar variations (quite apart from the familiar light-time effects).

TDB is defined in such a way that it keeps close to TT on the average, with the relativistic effects emerging as quasi-periodic differences of maximum amplitude rather less than 2 ms. This is negligible for many purposes, so that TT can act as a perfectly adequate surrogate for TDB in most cases, but unless taken into account would swamp long-term analysis of pulse arrival times from the millisecond pulsars.

Most of the variation between TDB and TT comes from the ellipticity of the Earth’s orbit; the TT clock’s speed and gravitational potential vary slightly during the course of the year, and as a consequence its rate as seen from an outside observer varies due to transverse Doppler effect and gravitational redshift. The main component is a sinusoidal variation of amplitude $0$$s$$.0017$; higher harmonics, and terms caused by Moon and planets, lie two orders of magnitude below this dominant annual term. Diurnal (topocentric) terms, a function of UT, are $2\mu$s or less.

The IAU 1976 resolution defined TDB by stipulating that TDB$-$TT consists of periodic terms only. This provided a good qualitative description, but turned out to contain hidden assumptions about the form of the solar-system ephemeris and hence lacked dynamical rigour. A later resolution, in 1991, introduced new coordinate time scales, TCG and TCB, and identified TDB as a linear transformation of one of them (TCB) with a rate chosen not to drift from TT on the average. Unfortunately even this improved definition has proved to contin ambiguities. The SLALIB sla_RCC function implements TDB in the way that is most consistent with the 1976 definition and with existing practice. It provides a model of TDB$-$TT accurate to a few nanoseconds.

Unlike TDB, the IAU 1991 coordinate time scales TCG and TCB (not supported by SLALIB functions at present) do not have their rates adjusted to track TT and consequently gain on TT and TDB, by about $0$$s$$.02$/year and $0$$s$$.5$/year respectively.

As already pointed out, the distinction between TT and TDB is of no practical importance for most purposes. For example when calling sla_PRENUT to generate a precession-nutation matrix, or when calling sla_EVP or sla_EPV to predict the Earth’s position and velocity, the time argument is strictly TDB, but TT is entirely adequate and will require much less computation.

The time scale used by the JPL solar-system ephemerides is called $T_{eph}$ and is numerically the same as TDB.

Predictions of topocentric solar-system phenomena such as occultations and eclipses require solar time UT as well as dynamical time. TT/TDB/ET is all that is required in order to compute the geocentric circumstances, but if horizon coordinates or geocentric parallax are to be tackled UT is also needed. A rough estimate of $\Delta T=ET-UT$ is available via the function sla_DT. For a given epoch (e.g. 1650) this returns an approximation to $\Delta T$ in seconds.

4.16 Calendars

The ordinary Gregorian Calendar Date, together with a time of day, can be used to express an epoch in any desired time scale. For many purposes, however, a continuous count of days is more convenient, and for this purpose the system of Julian Day Number can be used. JD zero is located about 7000 years ago, well before the historical era, and is formally defined in terms of Greenwich noon; for example Julian Day Number 2449444 began at noon on 1994 April 1. Julian Date is the same system but with a fractional part appended; Julian Date 2449443.5 was the midnight on which 1994 April 1 commenced. Because of the unwieldy size of Julian Dates and the awkwardness of the half-day offset, it is accepted practice to remove the leading ‘24’ and the trailing ‘.5’, producing what is called the Modified Julian Date: MJD = JD$-2400000.5$. SLALIB routines use MJD, as opposed to JD, throughout, largely to avoid loss of precision. 1994 April 1 commenced at MJD 49443.0.

Despite JD (and hence MJD) being defined in terms of (in effect) UT, the system can be used in conjunction with other time scales such as TAI, TT and TDB (and even sidereal time through the concept of Greenwich Sidereal Date). However, it is improper to express a UTC as a JD or MJD because of leap seconds.

SLALIB has six routines for converting to and from dates in the Gregorian calendar. The routines sla_CLDJ and sla_CALDJ both convert a calendar date into an MJD, the former interpreting years between 0 and 99 as 1st century and the latter as late 20th or early 21st century. The routines sla_DJCL and sla_DJCAL both convert an MJD into calendar year, month, day and fraction of a day; the latter performs rounding to a specified precision, important to avoid dates like ‘2005 04 01.***’ appearing in messages. Some of SLALIB’s low-precision ephemeris routines (sla_EARTH, sla_MOON and sla_ECOR) work in terms of year plus day-in-year (where day 1 = January 1st, at least for the modern era). This form of date can be generated by calling sla_CALYD (which defaults years 0-99 into 1950-2049) or sla_CLYD (which covers the full range from prehistoric times).

4.17 Geocentric Coordinates

The location of the observer on the Earth is significant in a number of ways. The most obvious, of course, is the effect of longitude and latitude on the observed $\left[Az,El\right]$ of a star. Less obvious is the need to allow for geocentric parallax when finding the Moon with a telescope (and when doing high-precision work involving the Sun or planets), and the need to correct observed radial velocities and apparent pulsar periods for the effects of the Earth’s rotation.

The SLALIB routine sla_OBS supplies details of groundbased observatories from an internal list. This is useful when writing applications that apply to more than one observatory; the user can enter a brief name, or browse through a list, and be spared the trouble of typing in the full latitude, longitude etc. The following Fortran code returns the full name, longitude and latitude of a specified observatory:

(Beware of the longitude sign convention, which is west +ve for historical reasons.) The following lists all the supported observatories:

The routine sla_GEOC converts a geodetic latitude (one referred to the local horizon) to a geocentric position, taking into account the Earth’s oblateness and also the height above sea level of the observer. The results are expressed in vector form, namely as the distance of the observer from the spin axis and equator respectively. The geocentric latitude can be found be evaluating ATAN2 of the two numbers. A full 3-D vector description of the position and velocity of the observer is available through the routine sla_PVOBS. For a specified geodetic latitude, height above sea level, and local sidereal time, sla_PVOBS generates a 6-element vector containing the position and velocity with respect to the true equator and equinox of date (i.e. compatible with apparent $\left[\alpha ,\delta \right]$). For some applications it will be necessary to convert to a mean $\left[\alpha ,\delta \right]$ frame (notably FK5, J2000) by multiplying elements 1-3 and 4-6 respectively with the appropriate precession matrix. (In theory an additional correction to the velocity vector is needed to allow for differential precession, but this correction is always negligible.)

See also the discussion of the routine sla_RVEROT, later.

4.18 Ephemerides

SLALIB includes routines for generating positions and velocities of Solar-System bodies. The accuracy objectives are modest, and the SLALIB facilities do not attempt to compete with precomputed ephemerides such as those provided by JPL, or with models containing thousands of terms. It is also worth noting that SLALIB’s very accurate star coordinate conversion routines are not strictly applicable to solar-system cases, though they are adequate for most practical purposes.

Earth/Sun ephemerides can be generated using the routines sla_EVP and sla_EPV, each of which predict Earth position and velocity with respect to both the solar-system barycentre and the Sun. The two routines offer different trade-offs between accuracy and execution time. For most purposes, sla_EVP is adequate: maximum velocity error is 0.42 metres per second; maximum heliocentric position error is 1600 km (equivalent to about $2$${}^{\prime }{}^{\prime }$ at 1 AU), with barycentric position errors about 4 times worse. The larger and slower sla_EPV delivers $3\sigma$ results of 0.005 metres per second in velocity and 15 km in position, and is particularly useful when predicting apparent directions of near-Earth objects. (The Sun’s position as seen from the Earth can, of course, be obtained simply by reversing the signs of the Cartesian components of the Earth : Sun vector.)

Geocentric Moon ephemerides are available from sla_DMOON, which predicts the Moon’s position and velocity with respect to the Earth’s centre. Direction accuracy is usually better than 10 km ($5$${}^{\prime }{}^{\prime }$) and distance accuracy a little worse.

Lower-precision but faster predictions for the Sun and Moon can be made by calling sla_EARTH and sla_MOON. Both are single precision and accept dates in the form of year, day-in-year and fraction of day (starting from a calendar date you need to call sla_CLYD or sla_CALYD to get the required year and day). The sla_EARTH routine returns the heliocentric position and velocity of the Earth’s centre for the mean equator and equinox of date. The accuracy is better than 20,000 km in position and 10 metres per second in speed. The position and velocity of the Moon with respect to the Earth’s centre for the mean equator and ecliptic of date can be obtained by calling sla_MOON. The positional accuracy is better than $30$${}^{\prime }{}^{\prime }$ in direction and 1000 km in distance.

Approximate ephemerides for all the major planets can be generated by calling sla_PLANET or sla_RDPLAN. These routines offer arcminute accuracy (much better for the inner planets and for Pluto) over a span of several millennia (but only $\pm 100$ years for Pluto). The routine sla_PLANET produces heliocentric position and velocity in the form of equatorial $\left[x,y,z,ẋ,ẏ,ż\right]$ for the mean equator and equinox of J2000. The vectors produced by sla_PLANET can be used in a variety of ways according to the requirements of the application concerned. The routine sla_RDPLAN uses sla_PLANET and sla_DMOON to deal with the common case of predicting a planet’s apparent $\left[\alpha ,\delta \right]$ and angular size as seen by a terrestrial observer.

Note that in predicting the position in the sky of a solar-system body it is necessary to allow for geocentric parallax. This correction is essential in the case of the Moon, where the observer’s position on the Earth can affect the Moon’s $\left[\alpha ,\delta \right]$ by up to $1^{\circ }$. The calculation can most conveniently be done by calling sla_PVOBS and subtracting the resulting 6-vector from the one produced by sla_DMOON, as is demonstrated by the following example:

The output produced is as follows:

(An easier but less instructive method of estimating the topocentric apparent place of the Moon is to call the routine sla_RDPLAN.)

As an example of using sla_PLANET, the following program estimates the geocentric separation between Venus and Jupiter during a close conjunction in 2 BC, which is a star-of-Bethlehem candidate:

The output produced (the Ephemeris Time on the day in question, and the closest approach in arcseconds) is as follows:

For comparison, accurate JPL predictions give a separation $8$${}^{\prime }{}^{\prime }$ less than the above estimate, occurring $30^m$ earlier (see Sky and Telescope, April 1987, p 357).

The following program demonstrates sla_RDPLAN.

Entering the following data (for 1927 June 29 at $5^{h}25^m$ ET and the position of Preston, UK):

produces the following report:

Inspection of the Sun and Moon data reveals that a total solar eclipse is in progress.

SLALIB also provides for the case where orbital elements (with respect to the J2000 equinox and ecliptic) are available. This allows predictions to be made for minor-planets and (if you ignore non-gravitational effects) comets. Furthermore, if major-planet elements for an epoch close to the date in question are available, more accurate predictions can be made than are offered by sla_RDPLAN and sla_PLANET.

The SLALIB planetary-prediction routines that work with orbital elements are sla_PLANTE (the orbital-elements equivalent of sla_RDPLAN), which predicts the topocentric $\left[\alpha ,\delta \right]$, and sla_PLANEL (the orbital-elements equivalent of sla_PLANET), which predicts the heliocentric $\left[x,y,z,ẋ,ẏ,ż\right]$ with respect to the J2000 equinox and equator. In addition, the routine sla_PV2EL does the inverse of sla_PLANEL, transforming $\left[x,y,z,ẋ,ẏ,ż\right]$ into osculating elements.

Osculating elements describe the unperturbed 2-body orbit. Depending on accuracy requirements, this unperturbed orbit is an adequate approximation to the actual orbit for a few weeks either side of the specified epoch, outside which perturbations due to the other bodies of the Solar System lead to increasing errors. Given a minor planet’s osculating elements for a particular date, predictions for a date only 100 days earlier or later are likely to be in error by several arcseconds. These errors can be reduced if new elements are generated which take account of the perturbations of the major planets, and this is what the routine sla_PERTEL does. Once sla_PERTEL has been called, to provide osculating elements close to the required date, the elements can be passed to sla_PLANEL or sla_PLANTE in the normal way. Predictions of arcsecond accuracy over a span of a decade or more are available using this technique.

Three different combinations of orbital elements are provided for, matching the usual conventions for major planets, minor planets and comets respectively. The choice is made through the argument JFORM:

The mean motion, $n$, tells sla_PLANEL the mass of the planet. If it is not available, it should be calculated from $n^2{a}^3=k^2\left(1+m\right)$, where $k=0.01720209895$ and m is the mass of the planet ($M_{\odot }=1$); $a$ is in AU.

Note that for any given problem there are up to three different epochs in play, and it is vital to distinguish clearly between them:

• The epoch of observation: the moment in time for which the position of the body is to be predicted.
• The epoch defining the position of the body: the moment in time at which, in the absence of purturbations, the specified position—mean longitude, mean anomaly, or perihelion—is reached.
• The epoch of osculation: the moment in time at which the given elements precisely specify the body’s position and velocity.

For the major-planet and minor-planet cases it is usual to make the epoch that defines the position of the body the same as the epoch of osculation. Thus, for planets (major and minor) only two different epochs are involved: the epoch of the elements and the epoch of observation. For comets, the epoch of perihelion fixes the position in the orbit and in general a different epoch of osculation will be chosen. Thus, for comets all three types of epoch are involved. How many of the three elements are present in a given SLALIB argument list depends on the routine concerned.

Two important sources for orbital elements are the Horizons service, operated by the Jet Propulsion Laboratory, Pasadena, and the Minor Planet Center, operated by the Center for Astrophysics, Harvard. The JPL elements (heliocentric, J2000 ecliptic and equinox) and MPC elements correspond to SLALIB arguments as shown in the following table, where “(rad)” means conversion from degrees to radians, and “(MJD)” means “subtract 2400000.5D0”:

Conventional elements are not the only way of specifying an orbit. The $\left[x,y,z,ẋ,ẏ,ż\right]$ state vector is an equally valid specification, and the so-called method of universal variables allows orbital calculations to be made directly, bypassing angular quantities and avoiding Kepler’s Equation. The universal-variables approach has various advantages, including better handling of near-parabolic cases and greater efficiency. SLALIB uses universal variables for its internal calculations and also offers a number of routines which applications can call.

The universal elements are the $\left[x,y,z,ẋ,ẏ,ż\right]$ and its epoch, plus the mass of the body. The SLALIB routines supplement these elements with certain redundant values in order to avoid unnecessary recomputation when the elements are next used.

The routines sla_EL2UE and sla_UE2EL transform conventional elements into the universal form and vice versa. The routine sla_PV2UE takes an $\left[x,y,z,ẋ,ẏ,ż\right]$ and forms the set of universal elements; sla_UE2PV takes a set of universal elements and predicts the $\left[x,y,z,ẋ,ẏ,ż\right]$ for a specified epoch. The routine sla_PERTUE provides updated universal elements, taking into account perturbations from the major planets. Starting with universal elements, the routine sla_PLANTU (the universal elements equivalent of sla_PLANTE) predicts topocentric $\left[\alpha ,\delta \right]$.

4.19 Radial Velocity and Light-Time Corrections

When publishing high-resolution spectral observations it is necessary to refer them to a specified standard of rest. This involves knowing the component in the direction of the source of the velocity of the observer. SLALIB provides a number of routines for this purpose, allowing observations to be referred to the Earth’s centre, the Sun, a Local Standard of Rest (either dynamical or kinematical), the centre of the Galaxy, and the mean motion of the Local Group.

The routine sla_RVEROT corrects for the diurnal rotation of the observer around the Earth’s axis. This is always less than 0.5 km/s.

No specific routine is provided to correct a radial velocity from geocentric to heliocentric, but this can easily be done by calling sla_EVP as follows (array declarations etc. omitted):

The maximum value of this correction is the Earth’s orbital speed of about 30 km/s. A related routine, sla_ECOR, computes the light-time correction with respect to the Sun. It would be used when reducing observations of a rapid variable-star for instance. For pulsar work the sla_EVP routine is not sufficiently accurate for phase predictions, being limited to about 25 ms. The alternative sla_EPV routine will deliver pulse arrival times accurate to 50 $\mu$s, but is significantly slower.

To remove the intrinsic $\sim 20$ km/s motion of the Sun relative to other stars in the solar neighbourhood, a velocity correction to a local standard of rest (LSR) is required. There are opportunities for mistakes here. There are two sorts of LSR, dynamical and kinematical, and multiple definitions exist for the latter. The dynamical LSR is a point near the Sun which is in a circular orbit around the Galactic centre; the Sun has a “peculiar” motion relative to the dynamical LSR. A kinematical LSR is the mean standard of rest of specified star catalogues or stellar populations, and its precise definition depends on which catalogues or populations were used and how the analysis was carried out. The Sun’s motion with respect to a kinematical LSR is called the “standard” solar motion. Radial velocity corrections to the dynamical LSR are produced by the routine sla_RVLSRD and to the adopted kinematical LSR by sla_RVLSRK. See the individual specifications for these routines for the precise definition of the LSR in each case.

For extragalactic sources, the centre of the Galaxy can be used as a standard of rest. The radial velocity correction from the dynamical LSR to the Galactic centre can be obtained by calling sla_RVGALC. Its maximum value is 220 km/s.

For very distant sources it is appropriate to work relative to the mean motion of the Local Group. The routine for computing the radial velocity correction in this case is sla_RVLG. Note that in this case the correction is with respect to the dynamical LSR, not the Galactic centre as might be expected. This conforms to the IAU definition, and confers immunity from revisions of the Galactic rotation speed.

4.20 Focal-Plane Astrometry

The relationship between the position of a star image in the focal plane of a telescope and the star’s celestial coordinates is usually described in terms of the tangent plane or gnomonic projection. This is the projection produced by a pin-hole camera and is a good approximation to the projection geometry of a traditional large f -ratio astrographic refractor. SLALIB includes a group of routines which transform star positions between their observed places on the celestial sphere and their $\left[x,y\right]$ coordinates in the tangent plane. The spherical coordinate system does not have to be $\left[\alpha ,\delta \right]$ but usually is. The so-called standard coordinates of a star are the tangent plane $\left[x,y\right]$, in radians, with respect to an origin at the tangent point, with the $y$-axis pointing north and the $x$-axis pointing east (in the direction of increasing $\alpha$). The factor relating the standard coordinates to the actual $\left[x,y\right]$ coordinates in, say, millimetres is simply the focal length of the telescope.

Given the $\left[\alpha ,\delta \right]$ of the plate centre (the tangent point) and the $\left[\alpha ,\delta \right]$ of a star within the field, the standard coordinates can be determined by calling sla_S2TP (single precision) or sla_DS2TP (double precision). The reverse transformation, where the $\left[x,y\right]$ is known and we wish to find the $\left[\alpha ,\delta \right]$, is carried out by calling sla_TP2S or sla_DTP2S. Occasionally we know the both the $\left[x,y\right]$ and the $\left[\alpha ,\delta \right]$ of a star and need to deduce the $\left[\alpha ,\delta \right]$ of the tangent point; this can be done by calling sla_TPS2C or sla_DTPS2C. (All of these transformations apply not just to $\left[\alpha ,\delta \right]$ but to other spherical coordinate systems, of course.) Equivalent (and faster) routines are provided which work directly in $\left[x,y,z\right]$ instead of spherical coordinates: sla_V2TP and sla_DV2TP, sla_TP2V and sla_DTP2V, sla_TPV2C and sla_DTPV2C.

Even at the best of times, the tangent plane projection is merely an approximation. Some telescopes and cameras exhibit considerable pincushion or barrel distortion and some have a curved focal surface. For example, neither Schmidt cameras nor (especially) large reflecting telescopes with wide-field corrector lenses are adequately modelled by tangent-plane geometry. In such cases, however, it is still possible to do most of the work using the (mathematically convenient) tangent-plane projection by inserting an extra step which applies or removes the distortion peculiar to the system concerned. A simple $r_1=r_0\left(1+K{r}_0^2\right)$ law works well in the majority of cases; $r_0$ is the radial distance in the tangent plane, $r_1$ is the radial distance after adding the distortion, and $K$ is a constant which depends on the telescope ($\theta$ is unaffected). The routine sla_PCD applies the distortion to an $\left[x,y\right]$ and sla_UNPCD removes it. For $\left[x,y\right]$ in radians, $K$ values range from $-1/3$ for the tiny amount of barrel distortion in Schmidt geometry to several hundred for the serious pincushion distortion produced by wide-field correctors in big reflecting telescopes (the AAT prime focus triplet corrector is about $K=+178.6$).

SLALIB includes a group of routines which can be put together to build a simple plate-reduction program. The heart of the group is sla_FITXY, which fits a linear model to relate two sets of $\left[x,y\right]$ coordinates, in the case of a plate reduction the measured positions of the images of a set of reference stars and the standard coordinates derived from their catalogue positions. The model is of the form:

where the p subscript indicates “predicted” coordinates (the model’s approximation to the ideal “expected” coordinates) and the m subscript indicates “measured coordinates”. The six coefficients a–f can optionally be constrained to represent a “solid body rotation” free of any squash or shear distortions. Without this constraint the model can, to some extent, accommodate effects like refraction, allowing mean places to be used directly and avoiding the extra complications of a full mean-apparent-observed transformation for each star. Having obtained the linear model, sla_PXY can be used to process the set of measured and expected coordinates, giving the predicted coordinates and determining the RMS residuals in x and y. The routine sla_XY2XY transforms one $\left[x,y\right]$ into another using the linear model. A model can be inverted by calling sla_INVF, and decomposed into zero points, scales, $x/y$ nonperpendicularity and orientation by calling sla_DCMPF.

4.21 Numerical Methods

SLALIB contains a small number of simple, general-purpose numerical-methods routines. They have no specific connection with positional astronomy but have proved useful in applications to do with simulation and fitting.

At the heart of many simulation programs is the generation of pseudo-random numbers, evenly distributed in a given range: sla_RANDOM does this. Pseudo-random normal deviates, or “Gaussian residuals”, are often required to simulate noise and can be generated by means of the function sla_GRESID. Neither routine will pass super-sophisticated statistical tests, but they work adequately for most practical purposes and avoid the need to call non-standard library routines peculiar to one sort of computer.

Applications which perform a least-squares fit using a traditional normal-equations methods can accomplish the required matrix-inversion by calling either sla_SMAT (single precision) or sla_DMAT (double). A generally better way to perform such fits is to use singular value decomposition. SLALIB provides a routine to do the decomposition itself, sla_SVD, and two routines to use the results: sla_SVDSOL generates the solution, and sla_SVDCOV produces the covariance matrix. A simple demonstration of the use of the SLALIB SVD routines is given below. It generates 500 simulated data points and fits them to a model which has 4 unknown coefficients. (The arrays in the example are sized to accept up to 1000 points and 20 unknowns.) The model is:

The test values for the four coefficients are $C_1=+50.0$, $C_2=-2.0$, $C_3=-10.0$ and $C_4=+25.0$. Gaussian noise, $\sigma =5.0$, is added to each “observation”.

The program produces output like the following:

In this above example, essentially identical results would be obtained if the more commonplace normal-equations method had been used, and the large $1000\times 20$ array would have been avoided. However, the SVD method comes into its own when the opportunity is taken to edit the W-matrix (the so-called “singular values”) in order to control possible ill-conditioning. The procedure involves replacing with zeroes any W-elements smaller than a nominated value, for example 0.001 times the largest W-element. Small W-elements indicate ill-conditioning, which in the case of the normal-equations method would produce spurious large coefficient values and possible arithmetic overflows. Using SVD, the effect on the solution of setting suspiciously small W-elements to zero is to restrain the offending coefficients from moving very far. The fact that action was taken can be reported to show the program user that something is amiss. Furthermore, if element W(J) was set to zero, the row numbers of the two biggest elements in the Jth column of the V-matrix identify the pair of solution coefficients that are dependent.

A more detailed description of SVD and its use in least-squares problems would be out of place here, and the reader is urged to refer to the relevant sections of the book Numerical Recipes (Press et al., Cambridge University Press, 1987).

The routines sla_COMBN and sla_PERMUT are useful for problems which involve combinations (different subsets) and permutations (different orders). Both return the next in a sequence of results, cycling through all the possible results as the routine is called repeatedly.