A&A 446, 747-771 (2006)
DOI: 10.1051/0004-6361:20053818
E. W. Greisen^{1} - M. R. Calabretta^{2} - F. G. Valdes^{3} - S. L. Allen^{4}
1 - National Radio Astronomy Observatory,
PO Box O,
Socorro, NM 87801-0387, USA
2 - Australia Telescope National Facility,
PO Box 76,
Epping, NSW 1710, Australia
3 - National Optical Astronomy Observatories,
PO Box 26732,
Tucson, AZ 85719, USA
4 - UCO/Lick Observatory,
University of California,
Santa Cruz, CA 95064, USA
Received 12 July 2005 / Accepted 5 October 2005
Abstract
Greisen & Calabretta (2002, A&A, 395, 1061) describe a generalized
method for specifying the coordinates of FITS data samples. Following
that general method, Calabretta & Greisen (2002, A&A, 395, 1077) describe
detailed conventions for defining celestial coordinates as they are
projected onto a two-dimensional plane. The present paper extends the
discussion to the spectral coordinates of wavelength, frequency, and
velocity. World coordinate functions are defined for spectral axes
sampled linearly in wavelength, frequency, or velocity, linearly in
the logarithm of wavelength or frequency, as projected by ideal
dispersing elements, and as specified by a lookup table.
Key words: methods: data analysis - techniques: image processing - techniques: radial velocities - techniques: spectroscopic - astronomical data bases: miscellaneous
The present paper is the third in a series of papers that establishes methods by which information about the physical coordinates of FITS data may be transferred along with the binary image, random groups, and table data. In "Paper I'' Greisen & Calabretta (2002) describe a revised method for transferring coordinate information in the FITS header and outline some rules governing the values assigned to the new standard header keywords. In "Paper II'' Calabretta & Greisen (2002) specify the conventions necessary to define celestial coordinates in a two-dimensional projection of the sky. This paper defines the parameters and conventions needed to specify spectral information including frequency, wavelength, and velocity. In addition to these basic conversions, a world coordinate function is defined to describe ideal optical dispersers of several common types. Several concepts that were suggested by spectroscopic issues are generalized to apply to all types of coordinate axes. These are a generalized description of non-linear algorithms, the -LOG and -TAB non-linear algorithms, and an axis-naming keyword. In "Paper IV'', Calabretta et al. (2005) specify methods to describe the distortions inherent in the image coordinate systems of real astronomical data.
Paper I describes the computation of the world or physical coordinates
as a multi-step process. The vector of pixel offsets from the
reference point is multiplied by a linear transformation matrix and
then scaled to physical units. Mathematically, this is given by
The final step in the computation is the conversion of these linear relative coordinates into the actual physical coordinates. The conventions to be applied to spectral axes, i.e. to frequency, wavelength, and velocity axes, are described in this paper.
Figure 1: Conventional velocities as a function of true radial velocity for selected values of the transverse velocity. The apparent radial velocity is plotted in the left panel, the radio velocity V/c is plotted in the center panel, and the optical velocity Z/c is plotted in the right panel. Note that, at a transverse velocity of 0.9 c, the red shift z always exceeds 1. Note that each family of curves intersects the axis of zero velocity measure at the same values of because the redshift is zero for these combinations of and . Thus, regardless of the value of , for any given redshift the velocity measures agree on whether the object appears to be receding or approaching. However, the object may appear to be receding at high speed even when it is actually approaching at high speed, though the converse is never true. At transverse velocities of all of the velocity measures are positive (receding) regardless of whether the object is actually receding or approaching. Furthermore, for below -0.5, the faster the object approaches, the smaller the transverse velocity required to make it appear to be receding! | |
Open with DEXTER |
At this stage in the discussion, we consider the spectral world coordinate at the reference point on all other axes. This covers the relatively simple case in which the spectral world coordinate is not dependent on the other world coordinates. Methods to describe deviations from this assumption due to the choice of spectral reference system are discussed in Sect. 7, while instrumental distortions are discussed briefly in Sect. 9 and, in more detail, in Paper IV.
Spectral coordinates are commonly given in units of frequency, wavelength, velocity, and other parameters proportional to these three. The coordinate types discussed here are then frequency, wavelength, and "apparent radial velocity'' denoted by the symbols , , and . There are also three conventional velocities frequently used in astronomy. These are the so-called "radio'' velocity, "optical'' velocity, and redshift, denoted here by V, Z, and z and given by , and z = Z / c. The velocities are defined so that an object receding from the observer has a positive velocity. We assume throughout that each image has at most one spectral coordinate axis; a similar limit on celestial coordinates was assumed in Paper II.
As discussed by Lindegren & Dravins (2003), the apparent radial velocity to be described here is itself a conventional velocity. The shift of a spectral line from its rest frequency is caused by a variety of effects, not just the Doppler shift due to a radial motion. The time dilation of an object moving with respect to the observer causes a spectral shift, even if that motion is entirely transverse. Gravitational fields at the radiating object and along the line of sight to the observer will shift the observed frequency. Furthermore, the apparent spectroscopic velocity may be produced at a peculiar point in the object, e.g. upwelling convective cells, rather than at the center of mass of the object. It may also be affected by absorption in intervening material and in particular by the cosmic redshift. We will use the term "apparent radial velocity'' here simply to refer to that pseudo-physical velocity described by the equations presented here. The apparent radial velocity then is a sum of all spectroscopic effects presented as if they were solely a radial velocity. While this may be sufficiently accurate for many uses, observers wishing to achieve very high accuracies should consult Lindegren & Dravins, and references therein, closely.
The effect of a transverse velocity deserves a little more discussion.
As given by Lang (1974),
For the purposes of this section and the next, we consider those spectral axes at a single celestial coordinate that are linearly sampled in wavelength, frequency, or apparent radial velocity. The radio and optical "velocities'' are directly proportional to frequency and wavelength, respectively, and thus do not constitute additional cases for the type of sampling. Frequency and wavelength axes may also be regularly sampled in their logarithm. Wavelengths are sometimes given "in air'' rather than in vacuum and denoted here by . This non-linear distinction is discussed in Sect. 4. Frequency may also be described in energy () units and in Kaysers ("wave number'' ) units.
Following the convention of Papers I and II, the first four characters of CTYPE specify^{} the coordinate type, the fifth character is `-' and the next three characters specify a predefined algorithm for computing the world coordinates from intermediate physical coordinates. When k is the spectral axis, the first four characters shall be one of the codes shown in Table 1. The mathematical symbols we will use for the codes and their default units are also listed. Units are specified with the CUNIT keyword. Standard values for unit strings are defined in Paper I. The IAU-standard prefixes for scaling the units are described in Paper I and should be used with all coordinate types, except that the dimensionless ones are not scaled.
Table 1: Spectral coordinate type codes, (characters 1-4 of CTYPE ).
The non-linear algorithm in use is specified by the final 3 characters of CTYPE . In spectral codes, the first of the three characters specifies the physical parameter type in which the data are regularly sampled, the second character is 2, and the third character specifies the physical parameter type in which the coordinate is expressed. Non-linear algorithm codes, the last four to be introduced later, are summarized in Table 2. When the algorithm in use is linear, the last four characters of CTYPE shall be blank. The original FITS paper by Wells et al. (1981) contained a number of suggestions for the values of CTYPE i. Those suggestions are to be superseded by the conventions of this paper and of Papers I and II.
Table 2: Non-linear algorithm codes, (characters 6-8 of CTYPE ).
The relationships between the basic physical quantities , and are shown in Table 3. The symbols and are the rest wavelength and frequency, respectively, of the spectral line used to associate velocity with observed wavelength and frequency. The relationships between the derived quantities and the basic quantities with which they are associated are shown in Table 4. Derivatives are included in both tables since they will be needed in the coordinate computation.
This section describes the computation of coordinates that are
sampled linearly or logarithmically in a coordinate of similar type,
and of coordinates that are sampled linearly in a coordinate of
different type. The intermediate world coordinate for spectral axis
k, given by Eq. (1), will be denoted, for convenience,
by
Linear coordinates are represented in CTYPE
by characters 1-4
containing any code in column one of Table 1 and
characters 5-8 blank. The world coordinate value for spectral axis
k is computed with the above definitions and
Eq. (1) as
Data are often sampled in logarithmic increments in one or more coordinates. For example, spectra are sometimes sampled in logarithmic increments of wavelength or frequency. With this type of sampling, the source motion, expressed as V, Z, or z, shifts the "pixel'' coordinates of spectral features by the same amount over the whole image. This allows relative velocities between spectra to be determined using cross-correlation methods in the pixel arrays.
For logarithmic axes, the last four characters shall be `-LOG'. While there are only three logarithmic coordinates
commonly used in spectroscopy, namely FREQ-LOG, WAVE-LOG
and AWAV-LOG, it would be unwise to forbid any coordinate
type with the -LOG non-linear algorithm. Many such combinations
may have little physical meaning or be intractable mathematically, but
these are simply reasons to be cautious when using -LOG. This
algorithm is evaluated simply
Table 3: Basic spectral transformation equations.
Table 4: Extended spectral transformation equations.
Logarithmic quantities are frequently expressed as log base 10 rather than as natural logarithms. For such data, the CDELT and CD k_ will need to compensate by including a factor of .
The generality of this algorithm, and the -TAB algorithm to
be introduced below, suggest the need for a more general description
of the coordinate than may be indicated in the first four letters in
the value of CTYPE
.
We hereby reserve the keyword
We now consider the case where an axis is linearly sampled in spectral variable X, but is to be expressed in terms of variable S.
Given the large number of spectral variables in Table 1 it is clear that there are very many pairwise combinations; in general, the relationship between any pair may be non-linear. However, each of the spectral variables in Table 1 is linearly related to one or other of , , or . Thus we may restrict X to one of these four basic variables.
Even with this restriction on X there are still many possible
combinations of X and S. In order to reduce the number still
further we introduce an intermediate variable, denoted by P,
that is the basic variable,
,
,
or
,
with which S is associated via a linear relation. This
associate variable is listed for each spectral variable in Col. 4 of
Table 1. Thus the sequence of transformations may be
summarized as
Table 3 lists the equations, X= X(P), and their inverses, P= P(X), linking the basic spectral types , and (discussion of is deferred to Sect. 4). These equations are generally non-linear. Likewise, Table 4 lists the linear relations, S= S(P), and their inverses, P= P(S), linking each spectral variable with its associate. When S is one of the basic types, , , or , then and S(P) is identity.
Thus the functional relationship between S and X is specified via intermediate variable P as S(X) = S(P(X)) with inverse X(S) = X(P(S)). Since S(P) is linear, Pand X must always differ, otherwise P(X) would be identity with the result that S(X) would be linear, contrary to the assumption of a non-linear axis.
Non-linear spectral coordinate codes are constructed on this basis by combining the spectral coordinate type code for S from Table 1 with the non-linear spectral algorithm code in Table 2. The first letter of the algorithm code defines X as frequency (F), wavelength (W), air wavelength (A), or apparent radial velocity (V), and the third letter likewise defines P. For example, in ZOPT-F2W, X is Frequency, P is Wavelength, and the non-linear conversion between the two (F2W) is defined in Table 3, Eq. (10). The desired spectral coordinate S is redshift (ZOPT), and this is related to associate variable P, the Wavelength, by the linear equation given in Table 4, Eq. (36).
It is apparent from the foregoing that it is possible to construct invalid spectral CTYPE coordinate codes. For example, ZOPT-F2V is unrecognized since z is associated with , not . That is not to say that this code could not be interpreted in principle - after all, equations linking z and could have been included in Table 4 - simply that it is not recognized in practice. The associate variable, P, was introduced in the first place to reduce the possible number of combinations and adding new ones like this would defeat that purpose. This is particularly relevant to the software implementation; ZOPT-F2V would tell software to chain its - function with a - z function but, in general, the latter function will not have been defined.
Consider now the computation of spectral coordinate S for an axis
that is linearly sampled in X. The statement that an axis is
linearly sampled in X simply means that
(43) |
(44) |
Equation (46) suggests a three-step algorithm chain:
(47) |
Table 5: Sample non-linear coordinate combinations, where w, the intermediate coordinate of Eq. (3), has a different meaning for each of these equations It has units given by the spectral coordinate type code as listed in Table 1.
As an example of a non-linear coordinate computation, consider a ZOPT-F2W axis. F in the F2W code indicates that the axis is linearly sampled in frequency, but it is to be expressed in terms of redshift as indicated by the spectral coordinate type of ZOPT. Table 1 indicates that redshift is associated with wavelength, hence the W in F2W is correct. Of the X, P, and S variables in the preceding section, the axis is linear in frequency so X is , the associated variable is wavelength, so P is , and the required variable is redshift, so S is z.
Since CRVAL
for the ZOPT-F2W axis would be recorded as a
redshift, ,
this must first be converted to frequency by
applying Eqs. (34),
One thing to notice about the equations of Table 5 is that the meaning of w differs for each. For example, Eq. (55) may not be obtained from Eq. (56) simply by multiplying both sides by c: in Eq. (56) w is a redshift (ZOPT), whereas in Eq. (55) it is an optical velocity (VOPT).
Aside from CRVAL
,
the coordinate computations of
Sect. 3.4 require one extra parameter when evaluating the
F2V, V2F, W2V, V2W, A2V, and
V2A non-linear algorithms, namely the rest frequency or
wavelength of the spectral-feature of interest. These are fundamental
physical parameters so, rather than use the PV i_ma parameters
defined in Paper I, which would tend to disguise them, the
special floating-valued keywords
Keyword RESTFREQ has been used in previous FITS files and should be recognized as equivalent to RESTFRQ.
The wavelengths so far discussed are measured in vacuum. However, dispersion coordinates for UV, optical, and IR spectra at nm are commonly given as wavelengths in air. The relative difference between the two varies between 0.028% and 0.032% across this range. To identify wavelengths measured in dry air at standard temperature and pressure rather than in vacuo, we introduce the coordinate type AWAV for which the reference value and increment must be expressed accordingly.
The two measures of wavelength are related by
Algorithm codes for the non-linear conversions are A2F, A2W, A2V, and their complements, F2A, W2A, and V2A as listed in Table 2. Use of the three-step procedure described in Sect. 3.4 would require as a function of each of the other spectral variables, together with their inverses and derivatives. Thus it is much simpler to handle air wavelengths as a separate, extra step in the algorithm chain. For example, to compute world coordinates for VRAD-A2V, the value of CRVAL would first have to be converted from radio velocity to air wavelengths via Eqs. (18), (10), and (67), and for Eq. (45) would be obtained from Eqs. (19), (7), and (66). Then, the value of computed for w via Eq. (42) would be transformed to vacuum wavelengths via Eq. (64), and then to radio velocity via Eqs. (6) and (20).
Figure 2: Geometry of gratings, prisms, and grisms. This simplified representation omits the collimation and focusing optics. Dashed lines mark ray paths in the plane of the figure - the "dispersion plane''. The normal to the grating/exit prism face and the normal to the detector plane are each projected onto the dispersion plane, and angles , , and are measured with respect to these projected normals. Usually the incident ray for a prism or grism is perpendicular to the entry face so that is equal to the prism angle, . Angle is wavelength-dependent, and consequently so is the offset in the dispersion direction on the detector. The intermediate spectral world coordinate, w, is proportional to . Reference wavelength follows the reference ray defined by and illuminates the reference point at . The normal to the detector plane is shown tilted by angle from the reference ray though typically this angle is zero. The grating spacing G^{-1} is indicated. | |
Open with DEXTER |
One common form of spectral data is produced by imaging the light from a disperser. The wavelengths of the light at some position in the image are related to the position and wavelength of the light in the field of view illuminating the disperser; the relation is defined by the physics of the disperser. In general the light received at a pixel in the detector is a superposition of different wavelengths from different points in the field of view. However, if the field of view is limited spatially, usually by aperture masks and fiber optics, each pixel receives light from only a small range of wavelengths. It is then possible to define a world coordinate function relating pixel position and effective wavelength. This is the basis of many astronomical spectrographs.
In the following section we use the physical relation applicable to the dispersers commonly used in astronomical spectrographs to define a world coordinate function for computing spectral coordinates. The relation applies to the simple, though common, case of single dispersers. More complex spectrographs with multiple dispersers, such as those using multiple passes through prisms, are not described by the methods of this section. The equations developed below also assume that the radiation enters perpendicular to the face of the prism, a condition not met by some widely used spectrometers. Alternatives to using this ideal world coordinate function, based on the physics of simple dispersers, are the table lookup described in Sect. 6 and empirical function fits provided by the the distortion function mechanism described in Paper IV.
We require that the dispersion occurs along just one direction on the detector and that the intermediate coordinates are computed so that only one world coordinate axis corresponds to wavelength. Paper IV describes how distortions and effects of the aperture shape can be removed to satisfy this requirement. The distortion correction is also used to remove aberrations causing departures from the ideal physical behavior of the dispersers assumed here.
The dispersers we consider are gratings, prisms, and grisms, which are a combination of a prism with a grating within or on the surface of the prism. Gratings may be reflecting or transmissive and may be fabricated with surface relief rulings, holographic surface relief patterns, and holographic volume phase patterns. In the discussion we use the terms grooves, lines, and ruling to refer to the periodic diffraction structures that produce the interference. By combining the two physical equations for interference and refraction a single world coordinate function covers all these cases with appropriate choice of parameters.
This section defines a world coordinate function using the basic laws of refraction and interference. It is beyond the scope of this discussion to give the full background for the laws and concepts underlying dispersive spectroscopy. Of the many texts on the subject, a standard reference is Astronomical Optics by Schroeder (1999).
Spectroscopic dispersers are based on interference and refraction effects which are naturally described in terms of wavelength. Moreover, they are most often used in the regime where wavelength is commonly the spectral coordinate. For these reasons the grism function is defined in terms of wavelength. The function assumes that all of the dispersion occurs at one place. This is why combinations of dispersers, other than a single grism, are not described by this function. This assumption means that for prisms and grisms the function is rigorously correct only when the light enters perpendicular to the face of the prism and the grating is at the exit face. However, even when these conditions are not exactly met the function may still be a good approximation with slight adjustments of the parameters.
The basic physical relationship between the wavelength, ,
the
angle of incidence of the light, ,
and the
diffracted/refracted angle, ,
is given by a combination
of the grating interference equation and Snell's law of refraction^{}:
By convention, the angles of incidence and diffraction/refraction are measured from the normal to the beam and is always measured positive in the anticlockwise direction. The sign changes on either side of the normal and this determines the sign of . Thus, for the reflection grating in Fig. 2, and , and vice versa for a transmission grating, prism, or grism.
Reflection grating geometry is sometimes defined by the angle measured from the camera axis to the collimator axis and the tilt tof the bisector relative to the grating normal. If and tobey the same sign convention as , then .
In spectrographs with prisms or grisms the requirement that the incident light be normal to the prism entry face means that is equal to the prism angle ( in Fig. 2). Even for oblique entry, identifying with the prism angle is often a good initial approximation.
The requirement that the diffraction and refraction occur at one point means that is fixed and independent of wavelength. The variation of with wavelength then defines the relation between wavelength and position on the detector. The reference wavelength is the wavelength at the reference point corresponding to the zero point of and the intermediate spectral world coordinate w.
The prism's dispersive power derives from the variation of its index
of refraction with wavelength. While this variation depends on the
material and is generally non-linear, it is sufficient to approximate
it by a first-order Taylor expansion about the reference wavelength,
:
In order to define a world coordinate function we need
as a
function of the intermediate world coordinate, w, which is
proportional to .
Since Eq. (70) gives
as a function of
it remains to determine the
relationship between
and .
First note that
Eq. (70) evaluated at
provides
Figure 2 shows angle ,
which is measured
from the reference ray to the camera axis using the same sign
convention as ,
i.e. if
is clockwise-positive as for a
grism then so is .
Normally
would be zero, but it is
included here to provide a small correction. Depending on the sign
convention for ,
simple geometry for a flat detector gives
In keeping with the preceding sections we wish to define general grating, prism, and grism world coordinate representations such as WAVE-GRI, FREQ-GRI, etc.
Bearing in mind that the grism equation, Eq. (70), is expressed in terms of wavelength, then given a FREQ-GRI axis, for example, it would be straightforward to convert the reference frequency, , given by CRVAL k, from frequency to wavelength via . However, it would not be valid to convert the intermediate spectral world coordinate, w, from frequency to wavelength like this because it is not a true frequency. While it may be a close approximation near the reference point it deviates at points away from it.
Thus, interpretation of an axis such as FREQ-GRI necessarily
involves a procedure similar to that of Sect. 3.4, and to
apply that methodology we need a parameter, the grism parameter,
,
that is a known function of the spectral variables and for
which the axis is linearly sampled; this will substitute for X in
Eq. (42). Since
is proportional to w
(73) |
Following Sect. 3.4 (with X replaced by ), we
have
As previously, we require that CDELT i or CD i_j be set so that
(79) |
(80) |
The foregoing provides everything needed to compute from w, and this forms the first step of a five-step algorithm chain for computing S from w:
(83) |
As usual, the inverse computation is effected by traversing the
algorithm chain in the reverse direction. The inverse equations
required for the backward steps have already been given except that
for Eq. (70):
Grism parameters required for these calculations are provided by the PV k_ma keywords defined in Table 6, with given by Eq. (71), and computed from the coordinate reference value, CRVAL k, by the appropriate equations from Tables 3 and 4. Default values for missing parameters are also defined in the table. Generally only the first three parameters will appear for gratings and the first five for grisms.
Table 6: Grism parameters, their default values, and required units.
A question raised above was how CDELT i or CD i_j may be set by a WCS composer so that .
w= s_{k} q_{k} | (85) |
It is the correct choice of sign for that resolves the sign ambiguity in Eq. (74).
Thus far we have ignored the distinction between vacuum and air wavelengths in the discussion of grism world coordinates. In fact, the variable that appears in the grism equation may be either, and in general in Eq. (68) is the index of refraction of the prism material with respect to the surrounding medium, either air or vacuum.
If the dispersion takes place in vacuum, as it does, for example, in a spectrograph on a spacecraft, then the grism equation is correct for vacuum wavelengths; if in air, it is correct for air wavelengths and in Sects. 5.1.1, 5.1.2, and 5.1.3 should be replaced everywhere by . Then as discussed in Sect. 4, the computation of quantities such as in Eq. (82) is best handled by using the vacuum wavelength, , as an intermediary, effectively introducing an extra step into the algorithm chain.
In order to distinguish between grisms operated in air and in vacuum we hereby reserve the use of GRI exclusively for vacuum operation and introduce GRA for operation in dry air at standard temperature and pressure.
Note that for real spectrographs operated in air at an observatory, the actual wavelength system that GRA describes is for air at the local conditions. The coordinate value and increment at the reference point are normally adjusted to the values at standard conditions during calibration. If the accuracy of the wavelength measurement requires it, any further correction between wavelength at local conditions and wavelength at standard conditions may be accomplished via the methods of Paper IV.
Figure 3: KPNO Coudé Feed spectrograph with a long focal length and a 3K CCD. | |
Open with DEXTER |
Figure 4: KPNO Hydra Fiber Bench Spectrograph using an echelle grating in order with a 2K CCD. | |
Open with DEXTER |
This section illustrates application of the GRA world coordinate function with three real-world examples of spectral images from Kitt Peak National Observatory (KPNO) spectrographs (Figs. 3-5). Each spectrum is of an arc calibration lamp that produces emission lines of known wavelength. The position of each line in the image along the spectral world coordinate axis is measured by centroiding on the spectral line profile and identified with the known rest wavelength to create a list of pixel positions and wavelengths.
Each example figure shows plots of the positions and wavelengths for a particular spectrograph. The pixel positions are plotted along the bottom axis. Rather than plot the wavelengths directly, where it would be difficult to see the shape of the curve, the difference or correction between the known (air) wavelengths and the simple linear world coordinate function AWAV ( ) are shown. The corrections are plotted along the left axes. In these examples the wavelength units are consistently Angstroms.
The top axis is labeled with simple AWAV linear coordinates. The right axis divides the wavelength correction by CDELT1, the linear dispersion at the reference pixel. This represents the shift, in pixels, of the wavelengths on the detector relative to where they would occur in a spectrograph with a linear dispersion relation.
Figure 5: KPNO MARS spectrograph with a 450 lines/mm volume phase holographic grism and a 2K CCD. | |
Open with DEXTER |
The departure of the data points from zero indicate the magnitude of the world coordinate errors that would occur by using the simple linear AWAV world coordinate function. The solid lines in the figures are the difference between the wavelengths produced by the AWAV-GRA grating coordinate function and the linear AWAV coordinate function evaluated at the pixel positions in the image.
The usefulness of the ideal grating coordinate function is that the curves go through the measured data points with an appropriate choice of parameters. The parameters are essentially those known for the spectrograph and disperser with small adjustments to some of the parameters to produce the best fit to the calibration data. (These small adjustments may be viewed as corrections for the simplifications in the optical model.) In these examples the dispersion function represented by the grating world coordinate relation is as good as typically provided by empirical polynomial functions. Higher order effects due to aberrations are corrected by the distortion correction methods defined in Paper IV.
A close reading of the equations above will reveal that the seven grism parameters listed in Table 6 are not independent. We have chosen these parameters because of their physical meaning. However, the independent parameters are and . It is these combinations of parameters which must be used in fitting data.
Each figure shows the arc line measurements, the coordinate function curve, and the relevant world coordinate keywords used to produce the curve. Not all of the WCS keywords are shown.
The first example shows the behavior of a long focal length spectrograph with a reflection grating used at a low angle of incidence. Because of the long focal length the deviation from linearity is relatively small though still clearly significant. The grating has a density of 316 lines/mm and is used in first order to produce a spectrum centered near 5225.2 Å with a dispersion of 0.43 Å/pixel.
The next example uses a 316 lines/mm echelle reflection grating, a grating designed for use at large angles from the grating normal, operated in a higher order. It is used to produce a spectrum centered near 5136.8 Å with a dispersion of 0.13 Å/pixel. The departure from linearity is in the opposite sense from the other examples because the echelle grating is used with the diffracted angle at a greater angle than the incidence angle. The higher order results in a fairly large departure from a linear WCS.
The last example illustrates use of a grism. However, this is not the more common grism with a ruled grating on the output face and the input face oriented normal to the incident beam. Instead, the grating is sandwiched in the middle of the prism. The prism is oriented with the beam entering and leaving the prism at equal angles to the faces resulting in a straight through configuration as with a standard grism.
Another unusual feature of this grism is that it uses a volume phase holographic (VPH) grating. While the intensity response is different from a surface relief grating (ruled or holographic) the dispersion behavior is the same.
The full optical equation is complex even in this symmetric configuration, but as shown in the figure using the GRA function, a very good description of the coordinates can be obtained.
The prism has a angle with an index of refraction of 1.764 near the reference wavelength. The grating has an equivalent interference pattern of 450 lines/mm. Using these values and adjusting the derivative of the index of refraction produces Fig. 5.
There are numerous instances in which a physical coordinate is well defined at each pixel along an image axis, but the relationship of the coordinate values between pixels cannot be described by a simple functional form. An obvious example of this would be a three-dimensional image consisting of a sequence of two-dimensional images of an astronomical object recorded at an arbitrary sequence of times determined in part by weather and time assignment committees. As another example, the calibration of some spectrographs, such as those employing a diode array detector, is represented best by a list of frequencies or wavelengths for each pixel on the spectral axis rather than some functional form.
To support such representations for primary images and image extensions, we define a table-lookup form for the value of CTYPE as xxxx-TAB, where xxxx are four letters representing the type of coordinate, e.g. TIME or FREQ. As in Paper I, which established the "4-3'' convention for CTYPE i, the coordinate xxxx is a "real'' coordinate, such as FREQ, not an intermediate coordinate, such as FREQ-F2W, requiring an additional linear or non-linear algorithm in order to be converted into a physical coordinate.
Consider first the case of a single (one-dimensional) coordinate axis. The -TAB algorithm uses a list of coordinate values, the coordinate array, to record coordinate values of the appropriate type for the coordinate axis. A second list of the same length, the indexing vector, may be used in addressing the coordinate array. The indexing vector provides one level of indirection which may be used to vary the sampling frequency of the coordinate array along the coordinate axis. The coordinate could then be sampled at smaller intervals over that portion of the range in which the instrument is more non-linear and sampled more coarsely over regions in which it is better behaved. The indexing algorithm, based on linear interpolation, is defined in Sect. 6.1.2. If the indexing vector is absent, it is taken to have values , where K is the number of values in the coordinate array. See Sect. 6.2.2 for an example of the use of the indexing vector.
The concept described above covers separable (one-dimensional) axes only. It may be extended simply to M non-separable axes so long as the indexing vectors for each of the M axes are separable. Mcoordinate values are required for each of the possible index positions. Therefore, the coordinates will be in a single (1+M)-dimensional array. This coordinate array will have dimensions , where K_{m} is the maximum value of the index on axis m+1 of the coordinate array. For simplicity, degenerate axes are forbidden; therefore, K_{m} > 1. The indexing vectors for each of the M axes each contain a one-dimensional array of length K_{m}.
Figure 6: -TAB logic flow with and without an indexing vector. The coordinate array subscript m associated with intermediate world coordinate axis i is specified with keyword PV i_. In the case of an independent -TAB axis it would have value 1. Note that is computed either with the CD i_j form of the linear transformation matrix or the PC i_j plus CDELT i form. | |
Open with DEXTER |
The data in the indexing vectors must be monotonically increasing or decreasing, although two adjacent index values in the vector may have the same value. However, it is not valid for an index value to appear more than twice in an index vector, nor for an index value to be repeated at the start or at the end of an index vector. Furthermore, repeated index values are allowed only in the case of one-dimensional separable axes. (See the following section for a discussion of the reasons for these limitations and how interpolation is done in such cases; see also Sect. 6.2.3 for an example of equal index values and further discussion of interpolation.) Application programs must not sort the index and coordinate arrays since this makes the relative order of the two equal index values indeterminate. The requirement for monotonic index values should eliminate any need for sorting. Note that it does not imply any ordering of the values in the coordinate array.
The indexing algorithm is illustrated schematically in
Fig. 6 for the one-dimensional case. In the general
case, to determine the M non-separable coordinate values C_{m}, one
first determines the M indices,
for addressing the
appropriate location in the table. If intermediate world coordinate
axis i is associated with the
axis in the coordinate
array, one evaluates
(87) |
In detail, the algorithm for computing
,
and thence
C_{m}, is as follows. Scan the indexing vector,
,
sequentially starting from the first element, ,
until
a successive pair of index values is found that encompass (i.e. such that
for monotonically
increasing index values or
for
monotonically decreasing index values for some k). Then, when
,
interpolate linearly on the indices
In the case where an index value is equal to , the algorithm above will find the interval for monotonically increasing index values or for monotonically decreasing index values, except when . Since two consecutive index values may be equal, the index following the matched index must be examined. If (or ), then and C_{m} are undefined.
Linear interpolation via Eq. (88) applies for in the range to inclusive. Outside this range, for K > 1, linear extrapolation is allowed for values of such that or for monotonic increasing index values, and for or for monotonic decreasing index values. Extrapolation is also allowed for K = 1with in the range (noting that should be equal to 1 in this case) whence .
The value of
derived from
must lie in the range
.
These extrapolation limits permit
assignment of coordinates to the regions of the pixels on the boundary
of the array which are outside of the centers of the boundary pixels
but within the conceptual "edge'' of the boundary pixels. In the
case of a single separable coordinate with
,
the coordinate value is given by
Conceptually, to compute the change in coordinate value between and , in the case when and/or , difference the two values of C_{m}obtained for and in the limit that goes to zero. In practice, the computation is straightforward and it is not necessary to take limits.
The inverse computation, in which one determines a given a coordinate C_{m}, is relatively straightforward, at least in the case of a single separable coordinate. One scans the coordinate vector from the start looking for a pair of values that encompass C_{m}. Having found a pair, one must then check the indexing vector. If the two index values are unequal, then the correct pair has been found. Otherwise the search must be continued. When a correct coordinate pair has been found, the inverse of Eq. (89) is applied to determine . The inverse of Eq. (88) then yields .
It is understood that the general keywords CUNIT , CRDER , and CSYER apply to the output coordinate C_{m} rather than the -TAB coordinate keywords CRVAL et al. Similarly, if the table lookup determines celestial coordinates, the general keywords RADESYS and EQUINOX apply to the output celestial coordinates rather than the input -TAB coordinates.
Standard FITS binary tables extensions (XTENSION = `BINTABLE') will be used to hold the coordinate array in a single cell of a column of a one-row table. The length of this array is given in the FITS table header by the repeat count in the TFORM n keyword, where n is the number of the column containing the coordinate array. The dimensions of the coordinate array will be given in the FITS table header by the keyword TDIM n set to `(M, K_{1}, K_{2}, K_{M})', where n is the column containing the coordinate array. Note in particular that if M = 1 the coordinate array is formally two-dimensional though the first axis is degenerate. The repeat count in the TFORM n keyword value is the product of M and all the K_{m}. The indexing vectors for each of the M axes, if present, will occupy separate columns, each containing a one-dimensional array of length K_{m}.
The BINTABLE extension containing the coordinate array must be in the same FITS file as the data that reference it.
The parameters required by -TAB are the table extension name (EXTNAME), the table version number (EXTVER), the table level number (EXTLEVEL), the column name for the coordinate array (TTYPEn_{1}), the column name for the indexing vector (TTYPEn_{2}), and the axis number m associated with intermediate world coordinate axis i in the coordinate array. The keywords used for this purpose are PS i_ , PV i_ , PV i_ , PS i_ , PS i_ , and PV i_ , respectively. These are summarized in Table 7. For images, PS i_ has no default; for tables a missing or blank PS i_ is taken to be the current table (see below). PV i_ and PV i_ both have a default value of 1. PS i_ never has a default; it must be present and be assigned a value actually occurring in table PS i_ . If PS i_ is missing or has a value of all blanks, the indexing vector is taken to be a list of integers from 1 to K_{m}and becomes a direct index of axis PV i_ in the array specified by PS i_ . If PV i_ is present and assigned a value, then that value must actually occur in table PS i_ . Note that the values given to PS i_ and PS i_ are not case sensitive since the FITS Standard (Hanisch et al. 2001) states that "String comparisons with the values of TTYPE n keywords should not be case sensitive.''
Table 7: Parameter keywords used for the -TAB algorithm.
The use of -TAB for M related axes requires the header to specify the array PS i_ to be the same for each of the M axes. The dimensions of this array must be given in the header as . These dimensions determine the maximum range of the array index (1 to K_{m}) for each axis. The indexing vectors for each of the M axes, if present, will occupy separate columns, each containing a one-dimensional array of length K_{m}. The M values of PV i_ must account for all M axes. If any of these conditions are not met, the result is undefined.
If a FITS file contains multiple XTENSION HDUs (header-data units) with the specified EXTNAME, EXTLEVEL, and EXTVER, then the result of the WCS table lookup is undefined. If the specified FITS BINTABLE contains no column, or multiple columns, with the specified TTYPE n, then the result of the WCS table lookup is undefined. The specified FITS BINTABLE must contain only one row.
No units conversions are to be performed. CUNIT i must be the same as TUNITn of the binary table, where n is the column number corresponding to PS i_ .
We recommend strongly that the value chosen for EXTNAME always begin with the four letters WCS-. If this is done, generic programs will recognize the table as part of the WCS portion of the data and will be less likely to separate the table from the rest of the WCS.
The -TAB implementation is more complicated than most other WCS conventions because the coordinate system is not completely defined by keywords in a single FITS header. Software that supports -TAB must be able to gather all the necessary WCS parameters that are in general distributed over two FITS HDUs and in the body of the WCS extension table. The producers of FITS data products should consider the capabilities of the likely recipients of their files when deciding whether or not to use the -TAB convention, and in general should use it only in cases where other simpler WCS conventions are not adequate.
Binary table extensions containing array data columns may need a table-lookup function for coordinate values. Seemingly the most convenient form of -TAB would be one in which the -TAB array(s) are also table cells in the same row as the data array. However, a separate coordinate table would be more economical if it applied to data arrays in multiple rows. The -TAB keywords for such tables are iCTYPna, iSn_0a, iVn_1a, iVn_2a, iSn_1a, and iSn_2a, where n is the column number of the array of data and i is the intermediate world coordinate axis number. If iSn_0a is missing or blank, the present binary table is taken to be the coordinate lookup table. In that case, the coordinate array(s) are taken to be single-cell arrays in the same row as the data array and keywords iVn_1a and iVn_2a are ignored.
Strictly speaking, the -TAB representation is not required for binary or ASCII table extensions containing only one data value per cell because each of the coordinate values associated with the datum may be stored in separate columns. However, -TAB does provide a convenient method of solving one of the problems pertaining to such tables, namely identifying the data column. If the table contains a column of data and another one or more columns of coordinate values pertaining to those data, then one can define the coordinates of the data column as being of type -TAB. In this case, the critical keywords are iCTYPna to declare the coordinate axis type and iSn_1a to identify the corresponding coordinate column. Since the data value column contains only one value, the coordinate column should only contain one value. Then, the usual coordinate keywords (jCRPXn / jCRPna, ijPCna, ijCDna, iCDLTn / iCDEna) may be omitted since their defaults yield the desired result that output index pixel equals the input data pixel (both of which are 1 in the assumed case). Note too that tables of this type - one value per table cell - may be in either ASCII or binary table form.
To illustrate the use of -TAB with two non-separable axes,
let
us consider a specific example of a four-dimensional image array.
Assume that axes 1 and 3 are handled either linearly or by one of the
non-linear single-axis cases (including -TAB). However,
assume that axes 2 and 4 require table lookup in a mutually dependent
fashion. Thus i = 2 for m = 1 and i = 4 for m = 2 and
the required keywords are
PS2_0 = `WCS-TAB' PS4_0 = `WCS-TAB' |
PV2_1 = 1 PV4_1 = 1 |
PV2_2 = 1 PV4_2 = 1 |
PS2_1 = `COORDS' PS4_1 = `COORDS' |
PS2_2 = `INDEX2' PS4_2 = `INDEX4' |
PV2_3 = 1 PV4_3 = 2 |
A real example that could use this particular algorithm is represented by a spectral image of a portion of the sky taken with a single radio telescope. The observer commands the telescope to point at a regular grid of coordinates, but, due to wind loading and other pointing errors, the telescope achieves the commanded positions only approximately. The actual celestial coordinates observed are, however, accurately measured. These data could then be represented usefully as a two-dimensional array of spectra, but accurate celestial coordinates for each spectrum could be found only by a 2-dimensional table lookup.
Since there has been no generally agreed upon FITS format for spectral data with explicit wavelengths assigned to each pixel, data providers have resorted to defining their own formats. Examples from the Hubble Space Telescope and the Very Large Array are shown below, recast into the -TAB algorithm.
Two types of Hubble Space Telescope data serve as examples of the -TAB algorithm. The two cases illustrate the evolution from simple images to the more powerful constructs provided by FITS extensions and binary tables. The purpose of discussing these formats is not to explain the HST formats, but to illustrate a couple of types of data that can be represented by the -TAB algorithm. Therefore, some details of these formats are ignored.
The early HST spectrographs, FOS and GHRS, adopted a format based only on simple FITS images. The basic concept is that the spectral flux values are given in one image and the vacuum wavelengths, in Angstroms, are given as the pixel values in another image. The two are associated by filenames. The filenames have the same basename, but different filename extensions for an exposure. To represent these spectra with a FITS WCS based on the -TAB method, while continuing to use an image representation for the spectral fluxes, one replaces the separate wavelength image with a table extension.
Table 8: WCS keywords in the spectral image for FOS and GHRS example.
Table 9: WCS keywords in the STIS table for two of the spectral columns.
Table 8 shows the minimum WCS keywords required in the primary image header. The PS1_0 keyword value defines the coordinate table to be in a binary table extension with name `WCS-TAB'. The coordinate table extension consists of one column, named `WAVELENGTH', and one row. The array values and length correspond to the original image format.
Other WCS keywords defined in Paper I would also be included as needed. The intermediate coordinate transformation must produce indexing values equal to the pixel coordinate in the image. Since the defaults for these keywords are defined in Paper I to produce pixel coordinates, it is not strictly necessary to include these keywords.
The second generation HST spectrograph, STIS, adopted a binary table format very close to one of those defined for the -TAB algorithm. In the STIS one-dimensional extracted echelle format each exposure is stored in a separate binary table. The extracted data from an exposure consist of a number of echelle orders each of which is stored in a separate row. The wavelengths and various associated "spectra'' are array elements in different columns. An associated spectrum is an array of a single type data quantity such as fluxes, errors, or data quality flags.
To convert this format to the -TAB representation only requires adding the 1CTYP n and 1S n_ keywords to specify the coordinate type and coordinate column name for each spectral column n. There are six types of spectral quantities for each order: gross, background, net, flux, error, and data quality flag. So there must be six pairs of keywords. Since the same wavelength column applies to each of the spectra, the keyword values are repeated six times, but with different column numbers in the keywords. Table 9 shows the WCS keywords for two of the spectral columns. The 1S n_0 keywords are omitted, signifying that the coordinate array column is in the same table as the spectra.
As noted for the FOS and GHRS format, additional WCS keywords are included as needed and the intermediate coordinate transformation values may be absent since this defines an indexing by pixel coordinate. In the STIS data, the units of the wavelength array are vacuum Angstroms so the 1CUNI n keyword would be required.
Table 10: Sample coordinate table. Each displayed column is actually a one-cell array of the (single-row) table.
Radio interferometry data are now best represented in binary tables with each row representing the visibility at one time on one antenna pair with a vector of data representing the complex values for all polarizations and frequencies observed. These telescopes often allow the user to observe regularly-sampled spectral channels starting from L arbitrarily chosen intermediate frequencies (IFs). Furthermore, this spectral data pattern may be used with one receiver and, a few minutes later, with other receivers at very different frequencies. Because of the great flexibility in the choice of the IFs and receivers, the frequencies present in the data may only be described by a table. Because that description is very repetitive, we choose to put it in a table separate from the visibility table.
We can construct a single column table of all frequencies observed. However, by using a second column for an indexing vector, we can take advantage of the regular sampling of the spectral channels. If we set iCRVLn, iCDLTn, and iCRPXn (for binary table format) or CRVAL i, CDELT i, and CRPIX i (for random groups and image formats) to one, then Table 10 may represent the actual frequencies, where is the frequency of channel 1 in IF and is the increment in frequency between the regularly sampled channels of IF . This example is shown graphically in Fig. 7.
In general, the FITS WCS papers do not consider the process by which the physical effects represented by the values in the FITS array have been measured. The measurable quantity may have been sampled over regions small enough to be considered as points in a widely separated grid of world coordinates. Alternatively, the measurable quantity may have been integrated across overlapping adjacent regions. The FITS WCS papers do not indicate whether it is valid to presume that the values in adjacent array elements represent physically adjacent quantities, nor do they indicate whether it is valid to interpolate the WCS across the extent of an array element.
Some world coordinate axes are defined to have integral values with conventional meanings (e.g. COMPLEX, STOKES). Interpolating a world coordinate value across the extent of a FITS array element along the direction of an axis which is conventionally defined to be integral would clearly be inappropriate. Nevertheless, for other types of world coordinate axes in traditional FITS arrays, the adjacent elements in the NAXIS-dimensional array can often be presumed to represent adjacent locations in a measurement continuum. This adjacency is relevant, for example, to image display programs when the image is magnified so that one FITS array element extends across multiple display pixels. In such cases a FITS display program may attempt to interpolate WCS values across the extent of an array element. With CTYPE k values other than -TAB this is commonly done by linear interpolation between the presumably adjacent array elements. However for coordinate axes where the CTYPE k is described by the xxxx-TAB non-linear algorithm there can be no presumption of adjacency.
Figure 7: Example taken from radio interferometry using -TAB with an indexing vector. The FITS keywords shown are suitable for the random groups format. The observation is made at a number of frequencies, with a number of regularly spaced spectral channels beginning from each of a number of arbitrary base frequencies . The use of an indexing vector reduces the number of values in the table from one array of 30 to two arrays of 10 each. In a real case the number of spectral channels would be significantly larger, making the space savings significant. In this example, pixel p_{j} = 6 produces . This lies at . The resulting coordinate is then which, as one would expect, equals . Note that this example involves only a single independent -TAB axis, so that PV i_ . | |
Open with DEXTER |
Table 11: Top: keywords in main FITS header for -TAB example, bottom: table keywords and sample arrays for that example.
Notwithstanding these caveats, the use of -TAB does permit specification of changes in the coordinate value across the extent of a single array element. Consider a sequence of two-dimensional images of the sky which have been re-gridded such that the pixel array for each image shares the same 2-D WCS for celestial coordinates. A single 3-D FITS array can hold such a sequence of 2-D images. If the same instrument produced each of the images on different dates then the third dimension might be purely temporal. Alternatively, if different instruments (radio, IR, optical, X-ray) produced each of the images then the third dimension might be spectral. In practice, however, images in such a spectral 3-D image will typically also have different observation dates, so the image will actually span four dimensions of world coordinates.
Akin to the example of a long slit spectrum in Paper II, the WCS of this four-dimensional case can be represented by the 3-D array. This fourth dimension could be represented with NAXIS = 4 and NAXIS4 = 1, or with NAXIS = 3 and WCSAXES = 4. Whether the sequence of images is temporal or spectral, it is certain that the exposure duration or the spectral bandpass have a non-zero extent and the -TAB representation does provide the ability to communicate the start and end times of an observation, or the minimum and maximum wavelengths of an observation. For this example the WCS keywords are shown in Table 11, where the keywords relating to axes 1 and 2 are omitted because they are simply celestial coordinates as described by Paper II. In this example, note the clever use of the PC matrix to cause the fourth coordinate axis to depend on the third pixel axis, but not the degenerate fourth pixel axis. FITS writers can use details such as these to communicate how a display program might provide meaningful world coordinates at locations across an array element^{}.
In order to provide coordinate values across the entire array element, the coordinate values on the boundaries between the array elements are multiply defined. Such a result should be expected for array elements representing non-adjacent physical values. It is evident in this case that the order of the index arrays should not be permuted during the coordinate lookup. The lower portion of Table 11 shows a subset of the content of the table(s) referred to by PS _ and PS _. The four columns in the figure are meant to illustrate the association of table header keywords and table values. In the FITS file, the keywords would appear in the table header along with many other keywords while the four arrays of values would appear in four cells of the one-row binary table.
The indexing vectors in the table of Table 11
illustrate the reason to allow two values within the indexing
vector to be identical. The rules for the linear interpolation within
such vectors are clear. One must select the two index locations
having values immediately surrounding .
Thus in the present
example, if the time coordinate has
,
the third and
fourth index and coordinate values are interpolated. The output time
coordinate is
Frequencies, wavelengths, and apparent radial velocities are always referred to some selected standard of rest (reference frame), and while they are measured, of necessity, in the observer's rest frame, they may also be corrected to some other frame^{}. The velocity correction is computed from the vector dot product of the direction vector and the relative velocity vector of the two reference frames. In addition, the corrections from topocentric, the frame in which the observations are usually made, to geocentric and then to barycentric or heliocentric are dependent on the dot product with time-variable velocity vectors. As a consequence, the "velocity'' with respect to the reference frame depends on direction. Differential effects over a field of view may be important in some high precision applications; for example, they may amount to 5 km s^{-1} for the Local Group correction over a field of one degree. In another example, observations of Galactic CO over two-degree fields separated by two degrees failed to align by two spectral channels (Mangum et al. 2000).
Table 12: Recognized values for SPECSYS a, SSYSOBS a, and SSYSSRC a.
If a frequency, wavelength or apparent radial velocity associated with an image plane has been corrected to some other frame by transforming the CRVAL et al. values, then differential errors may arise at points away from the reference point. For example, it is normal in radio astronomy to observe a field while holding constant the velocity of the reference point with respect to the local standard of rest. See Sect. 10.2 for a more detailed discussion of the complications inherent in such data.
Nonetheless, each image plane shares a constant topocentric frequency
(or apparent radial velocity). The velocity (and
frequency and wavelength) with respect to the local standard of rest
is then a function of celestial coordinate within the image. In order
to denote this, we introduce two character-valued keywords. The
first,
Many spectrometers suffer from a significant variation in their spectral coordinate as a function of celestial coordinate. Such variations are the primary subject of Paper IV, while the concepts of this work refer to images produced by, or corrected to, an ideal spectrometer. In general we recommend that the frame corrections will get folded into the instrumental ones and that the image produced after the Paper IV corrections will have identical SSYSOBS and SPECSYS. This is not required, but would certainly simplify matters.
It is not appropriate for this work to define the parameters of
the recognized rest frames. The parameters needed to compute
geocentric frequencies/velocities from topocentric are the sidereal
time (or Earth rotation angle) and the observatory location. The
observing date is needed to convert from geocentric to barycentric
coordinates. The new keywords
Astronomers have traditionally quoted telescope positions using terrestrial longitude, latitude, and height, but, for computation of topocentric velocity, these values are immediately converted into a geocentric Cartesian triple. A precise conversion from longitude, latitude, and height, however, requires knowledge of the geodetic datum in which they are expressed. There are about 1000 geodetic datums in use. These are described by the United States Department of Defense (2000) in a printed document also available over the Internet. Some of the datums deviate from geocentricity by over a kilometer.
Table 13: New spectral keywords including forms for use in tables. All table keywords except the coordinate name have the same form for the BINTABLE vector image format and the pixel list format.
Rather than burden FITS with these details of geodetic tradition,
three new keywords
It may be helpful for the FITS writer to provide information on the
relative radial velocity between the observer and the selected
standard of rest in the direction of the celestial reference
coordinate. The keyword
The frame of rest defined with respect to the source (SOURCE) is
useful for comparing internal apparent radial velocities in
objects at different redshifts. This allows the FITS writer to apply
all the cosmological and other model-dependent corrections,
leaving a coordinate local to the object. For this frame of
rest, it is necessary to define the velocity with respect to some
other frame of rest. The keywords
The use of coordinate keywords to describe a multi-dimensional vector in a single element of a FITS binary table (Cotton et al. 1995), and a tabulated list of pixel coordinates in a FITS ASCII (Harten et al. 1988) or binary table was discussed in Papers I and II. In this section, we extend those discussions to the keywords specific to spectral coordinates. This convention is considered an integral part of the full world coordinates convention.
Table 13 lists the corresponding set of coordinate system keywords for use with each type of FITS image representation. The keywords defined in this paper and their allowed values are applicable to all three image storage formats. The following notes apply to the naming conventions used in Table 13:
There are instruments having spectral coordinates that are a function of other coordinates in the resulting image, e.g. celestial position. For example, optical astronomers use imaging Michelson interferometers and scanning Fabry-Perot instruments to produce three-dimensional images with two celestial and one spectroscopic coordinate. Because the path length through such instruments increases off the optical axis of the telescope, the observed wavelength is a function of celestial coordinate in each plane of the image cube. In general, this is a subject reserved for Paper IV since such instruments almost always have lesser distortions that need to be parameterized and corrected along with the path length correction.
One particularly difficult case should be mentioned here. Objective prism instruments produce data in which the image of each object in the field has been spread out into a spectrum. Thus the output image is the convolution of the sky with both a spatial point-spread-function and a spectral dispersive point-spread function, both in the two dimensions of the output image. The methods used to specify coordinates described in the present paper and in Papers I and II are unable to handle such complex data. However, rather than ignore the world coordinates entirely, writers of such images should consider providing the celestial coordinates that would be correct for the image if it had been taken at a single representative wavelength. That wavelength would be specified in a third WCS axis. Alternate axis descriptions could then be used to specify the same information for other wavelengths. Multi-order echelle spectrographs also produce images in which multiple spectral orders overlap each other as well as the celestial axes and, hence, their coordinates cannot be described by the methods of Papers II and III.
Analysis of both objective prism and echelle spectrograph images can produce tables of spectra that are easily described by the present methods.
Table 14: Example image header.
Table 15: Additional keywords for example header.
Let us consider an example. The partial FITS header in Table 14 was produced for data from the Very Large Array. The FITS WCS standard states that all WCS keywords must be reproduced for the alternate descriptions, but those for non-spectral axes have been omitted for the sake of brevity and clarity. We have also updated the keywords to the standards of Papers I, II, and III. We will derive additional header keywords that could also have been written. These are shown in Table 15, again leaving out axes 1 and 2 of the alternate WCS versions.
The spectral-line cube is regularly sampled in frequency and the primary description is of a linear frequency axis (from the CTYPE3 keyword value `FREQ'). CRVAL3 and CRPIX3 record that the 63-channel spectrum was centered with channel 32 at 1.37835117405 GHz, some 42 MHz from the rest frequency (RESTFRQ = 1.420405752 GHz) of neutral hydrogen (H I), the line under investigation. As indicated by the SPECSYS keyword, it is shown as if it had been observed at a constant topocentric frequency. For the moment, let us assume that the observations were made over a short interval and that this is correct.
As allowed by Paper I, we use alternate axis description letter
codes in a mnemonic fashion. Using Z for optical, the first
alternate axis description for the spectral axis reflects the
observer's request that a Doppler shift be applied to the H I line
frequency appropriate for a source with a barycentric
optical-convention velocity of 9120 km s^{-1}. Combining
Eqs. (30) and (6), we get
With CRVALZ set to
m s^{-1} it remains to compute
CDELT3Z. This is obtained by transforming CDELT3 from the
topocentric to the barycentric frame with Eq. (8) expressed
as
The correction for the observatory's motion due to the Earth's rotation and orbital motion with respect to the barycenter amounted to 26.108 km s^{-1} in the direction of the source. As illustrated in Table 15, this velocity shift can be represented in the FITS header with the new keyword VELOSYS.
The barycentric frequency may be used as another alternate axis description. The reference frequency and frequency increment for SPECSYSF = `BARYCENT' have already been computed above. The relevant keywords have axis description F (for frequency) in Table 15.
The frequency description may be translated into a wavelength description simply using Eqs. (10) and (11), shown in Table 15 using axis description W (for wavelength).
The frequency description with respect to the BARYCENT may also be expressed simply as a radio velocity (also regularly sampled) using Eqs. (20) and (21). This is shown as axis description version R (for radio).
A apparent radial velocity description (version V for velocity) requires the use of Eq. (14) for the reference value and Eq. (15) for the increment.
An assumption made throughout this paper, particularly in Sect. 7, is that data are provided with at least celestial coordinate information in addition to the spectral coordinates discussed here. The normal method of providing such data is with the full set of coordinate keywords including WCSAXES with a value of at least 3 (for two celestial and one spectral axis) even if only one celestial position was observed. For spectra in tables, it would be possible to have simple columns for the celestial coordinates, labeled as such. But, as discussed at the end of Sect. 6.1.4, this simple case may still use the full WCS nomenclature in order to associate the coordinate columns with the data column.
One important assumption made in Sect. 10.1 was that the topocentric velocity does not change significantly during the course of the observation. We now consider the consequences if this assumption is violated. In this we are interested in particular with three-dimensional "data cubes'' containing multiple samples on each of two spatial axes and one spectral axis.
For some instruments, it may be possible to correct the data for each short exposure to a consistent WCS and then combine exposures to improve the sensitivity. For example, single-dish radio telescopes observe a single direction at a time, allowing the spectral axis in each observed spectrum to be scaled and resampled before it is added into a three-dimensional image. For this type of instrument it is a relatively simple matter to correct for topocentric motion.
However, other instruments observe multiple directions simultaneously and also routinely combine multiple exposures to form an image. While the error over a single exposure may be negligible, images made by aperture synthesis radio telescopes, for example, are often based on multiple exposures taken months apart. These exposures are usually made with different array configurations, which may change on a timescale of months, in order to sample the Fourier space adequately for imaging and deconvolution. Therefore, these data are affected not only by the rotation of the Earth but also by its motion around the Sun which usually is not negligible.
As a specific example, the actual observation on which Sect. 10.1 is based required an integration time of several hours. During this time the observing frequency was adjusted for the topocentric Doppler shift so that the barycentric frequency at the map center remained constant. This ensured that the spectral feature of interest, associated with a source near the map center, remained in the same spectral channel. However, this means that the primary spectral axis description, based as it is in the topocentric frame, is not strictly correct. It was used in this data as an aid to interferometric analysis, the resulting error being considered negligible, but strictly speaking a barycentric frequency reference frame, the F alternative WCS in the header, is the correct one.
Even so, the barycentric frame is only strictly correct at the map center since the topocentric Doppler correction varies across the field-of-view; the correction is computed as the dot product of the observatory's velocity vector and a unit direction vector which, of course, varies across the field. This gives rise to a differential error that is greatest when the velocity vector is orthogonal to the direction to the source and least when it is parallel or anti-parallel. In other words, the differential correction is greatest when the correction is zero, and least when the correction is greatest.
For diurnal rotation the worst case differential error occurs for a telescope at the equator observing a source on the local meridian. For an angular offset, , from the tracking center the maximum differential error is which amounts to 8.7 m s^{-1}for . The worst case change in the differential error over the course of 12 hours occurs for an equatorial telescope observing a field at the north or south point on the horizon; it is or 17.5 m s^{-1} for . If data from exposures separated by several months are combined then the Earth's 30 km s^{-1} orbital velocity scales up the worst-case error to 1.0 km s^{-1} for . Note that, in the example of Sect. 10.1, the barycentric correction is close to the maximum value so the differential correction is nearly at its minimum.
The relevance of the SSYSOBSa keyword immediately becomes apparent in this context. It records the reference frame in which the frequency has no spatial variation; in other frames the reference frequency and frequency increment are only strictly correct at the reference point. While SSYSOBSa will most often indicate the topocentric frame, it is possible to regrid (resample) the spectral data by applying a position-dependent shift-and-scale so that SSYSOBSa has some other value.
Differential Doppler effects are often neglected for contemporary observations made with small angular fields or modest spectral resolution. However, instruments now under development will be capable of significantly greater sensitivity over wider fields-of-view and/or spectral resolution and hence will need to consider the effects discussed here.
The new keywords required for spectral coordinate systems are summarized in Table 13. SPECSYSa, SSYSOBSa, SSYSSRCa, and ZSOURCEa are used to specify velocity reference frames; MJD-AVG, DATE-AVG, OBSGEO-X, OBSGEO-Y, OBSGEO-Z, VELOSYSa, and ZSOURCEa enable conversion between these standards of rest; and RESTFRQa, RESTFREQ, and RESTWAVa define the spectral line for which velocities are measured. Variants of these keywords for use with tabular data are also defined.
These new keywords and their allowed values, along with new values for the standard header keywords formalized by Paper I, and the associated algorithms and methods introduced here allow an accurate description of spectral coordinates in FITS images. Wavelength, frequency, apparent radial velocity and oft-used quantities that are linearly related to them, such as redshift, energy, and radio velocity, may now be expressed as functions of pixel coordinate along an axis regularly sampled in wavelength, frequency, or apparent radial velocity.
A non-linear algorithm is also provided for spectral dispersers commonly used at optical wavelengths: gratings, prisms, and the combination of the two, grisms. Although developed for ideal dispersers, it was shown to apply quite well to real-world dispersers with suitable fine adjustment of the instrumental parameters.
The multi-dimensional -TAB table lookup algorithm developed for wavelength calibration will be useful for much more than just spectral axes. Several others of the methods introduced here are also completely general; logarithmic (-LOG) coordinates, and the CNAME keyword. The mathematical formalism described in Sect. 3.4 for constructing one-dimensional non-linear coordinate systems also has wide validity.
The discussion of Sect. 2, particularly Fig. 1, indicates the importance of transverse, as well as radial velocities, in spectroscopy at significant velocities. To be concrete, consider a relativistic jet emanating from a distant galaxy. In principle, the galaxy's systemic, cosmological redshift can be measured separately and used to correct the jet's observed redshift thereby providing its kinematic redshift (i.e. associated with a true velocity) in the reference frame of the galaxy.
Clearly knowledge of the velocity is fundamental in studying jet kinematics and dynamics. However, it can only be computed from the kinematic redshift if the jet's orientation angle is known. Note that the equations involving velocity in Table 3 are actually only correct if the transverse velocity is zero. For example, Eq. (12) is just a special case of Eq. (2). This is why the velocity is always referred to in the paper as being "apparent.''
There are instances where the orientation angle may be inferred by
geometry (e.g. by the observed tilt of an accretion disk) or by
modelling. Defining the orientation angle as
we have
(A.1) | |||
(A.2) | |||
(A.3) |
Table A.1: Velocity equations using orientation and total velocity.
There are other combinations that are of interest, including
It has been proposed that we introduce a single keyword to express
so that we may express velocities in terms of a true space
velocity
rather than an apparent velocity
as done
throughout this paper. However, astrophysical jets are thought to be
conical in form, to interact along their boundaries with the external
medium, and to change direction, even assuming helical shapes. Thus,
there is no single angle
that can normally be used to
describe all directions in an image and one should probably not even
assume that
is independent of redshift in any one direction.
For those cases in which a single angle may be appropriate, such as a
spectrum of a spatially limited region in a relativistic jet, we
reserve the keyword
Acknowledgements
The authors wish to acknowledge constructive comments and contributions from the following: Wim Brouw (Australia Telescope National Facility), Bob Hanisch (Space Telescope Science Institute), Jonathan McDowell (Harvard University), William Pence (Goddard Space Flight Center), Arnold Rots (CfA Harvard University), William Thompson (Goddard Space Flight Center), Doug Tody (National Radio Astronomy Observatory), and Patrick Wallace (UK Starlink).
The SLALIB software library developed by Patrick Wallace was helpful in the preparation of Sect. 7 and should be of considerable help to users of conventions described in this work. It is available via http://www.starlink.ac.uk/.
The National Radio Astronomy Observatory is a facility of the (U.S.) National Science Foundation operated under cooperative agreement by Associated Universities, Inc.
The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO.
The National Optical Astronomy Observatory is a facility of the (U.S.) National Science Foundation operated under cooperative agreement by Associated Universities for Research in Astronomy, Inc.
UCO/Lick Observatory is operated by the University of California.