Figure 9: Spectral data layout.
Flat-field division should be relatively straightforward; you take your raw data and divide it pixel for pixel by your flat-field. In practice, this may be acceptable for image data, but there are additional considerations to be taken into account when the data involved are spectral.
In principle, you can simply divide your image by the appropriate flat-field and use the result. For example:
In practice, you really don’t usually want to divide your data by the (typically) very large numbers found in flat-fields. Ideally, your flat-field should first be reduced to numbers around unity. The easiest way to do that is to divide the flat-field by its mean value, and the easiest way to find the mean value is with ‘istat’. The command
will give you the mean value over the whole image. If your flat-field has strange effects round the edges, you may prefer to limit the range of x and y values used by ‘istat’. For example, if your image is 800 by 800,
will only look at the central 600 by 600 part of the image. ‘istat’ prints out the value of the mean, and you can then divide the flat-field by that. For example, supposing the mean value were 350, the command
will generate a new flat-field (ffdiv) that will have a mean value around unity, as desired. You may find it convenient to know that ‘istat’ sets its output parameter ‘stat_mean’ to the mean value. So you could use
Figaro expects you to have your spectral data oriented so that the x dimension is the wavelength dimension. If this is not the case, use ‘irot90’ to switch the axes. In theory, Figaro does not care whether wavelength increases or decreases with x value, but in practice routines tend to be tested with data whose wavelength increases with x value, and odd bugs may turn up with ‘reversed’ data. You are recommended to get your data into the more common form, if necessary using ‘irevx’. Note that, as always in Figaro, ‘x’ means the horizontal axis when data is displayed using ‘image’.
The main problem with spectral flat-fields is that the flat-field data will usually vary in x, not because of the instrumental response, but simply because of the spectral response of the flat-field lamp. The usual way of dealing with this is to fit this spectral response—obtained by collapsing the flat-field in y—and then multiplying by the fitted value before dividing by the flat-field. Another way to look at this is to consider the result of dividing the flat-field by the fit to the spectral response. The result would be an image which was essentially flat in X, but which is ‘humped’ in Y—since most flat-fields fall off at the edges in Y due to instrumental vignetting. Dividing by this—in a two-dimensional manner—will give you images where the pixel to pixel variations in the detector have been corrected, along with any spatial vignetting of the instrument, but where the overall wavelength response of the instrument is not corrected.
(Of course, if you think your flat-field lamp is flat—which may be true at high dispersion—you can just divide by the raw flat-field and you’re taking out the instrumental wavelength response as well. However, the better way to deal with that is with your standard stars.)
The only problem, usually, is deciding on what represents an acceptable fit to the collapsed flat-field. One approach is to fit a polynomial to it. This at least has the advantage of being automatic, and if the result is satisfactory this is probably the best thing to do.
So this recipe is as follows:
The individual Figaro commands to do this might be the following
which, given a flat-field exposure in ‘ff’, and an image in ‘mydata’, performs a flat-field calibration as described above and puts the result in ‘mydataff’.
This is a sufficiently common way of proceeding that there is a single Figaro command that performs all of this, without the overheads of running four programs and generating the intermediate files. This is the ‘ff’ command. The command
is functionally equivalent to the above sequence. However, the advantage of splitting up the process is that you can compare ‘spres’ and ‘spfit’ to see how good the polynomial fit to the collapsed data really is. E.g. by
(You also have the advantage that you can control the limits used when collapsing the flat-field; you do not have to use ‘min’ and ‘max’ as the limits, although that is what ‘ff’ does.)
If the fit is not satisfactory, what can you do instead? Well, the problem is essentially that of generating a smoothed spectrum, and Figaro has a number of ways of doing that.
One of these should enable you to get a satisfactory result, although they all require more effort than does the simple ‘ff’.
The atmospheric features that mess up the red end of the spectrum may be calibrated out by multiplying by a calibration spectrum obtained from an observation of an object that is featureless in this region. A B star is usually used, hence the term ‘B-star calibration’.
The essential Figaro command for this calibration is ‘bsmult’. This multiplies the spectrum to be corrected by a B-star calibration spectrum—which is essentially obtained by fitting the continuum of a B star and then dividing that fitted continuum by the B-star observation. ‘bsmult’ allows for the difference in air masses of the two observations, which is what makes it necessary, rather than, say, the simpler ‘imult’.
‘bsmult’ is a fairly straightforward operation. The problems come in generating the calibration spectrum. At present, the simplest way is to display the B-star spectrum using ‘splot’ and then generate the continuum by hand using ‘cfit’. An alternative is to generate a mask for the lines in the spectrum using ‘mask’ and then fit the masked continuum using ‘mcfit’. The main problem with the latter approach is that the chances are that very little of the spectrum is actually uncontaminated either by stellar lines or by the atmospheric bands.
(The best solution is use an automatic program that fits splines between points on the spectrum that are known to be uncontaminated, but such a program is not available yet—nor is a really good list of such points that will apply for all wavelength ranges and dispersions.)
The calibration spectrum should really be exactly 1.0 at all points not affected by the atmospheric bands, and it is probably worth displaying the calibration spectrum using ‘splot’ and then using ‘cset’ to set such regions to 1.0 interactively. (This is something else that could be made automatic eventually. Fortunately, it isn’t necessary to do this very often.)
Note that ‘bsmult’ requires that both the calibration and the object being corrected have valid air masses. Air masses are stored in NDF data files as ‘file.MORE.FIGARO.SECZ’, and in DST data files as ‘file.OBS.SECZ’. If necessary they may be set by hand using, e.g.
Also note that the correction applied by ‘bsmult’ is multiplicative—this means that it is not suitable for data that is in logarithmic flux units, such as magnitudes.
If the response of a filter has been tabulated, then the table of values may be used to generate a ‘spiketrum’, which may then be turned into a calibration spectrum. See Section 4.4 for details about these things, which are essentially ways of turning sets of tabulated values into spectra by interpolation.
Say the spectrum to be corrected for a filter response is called ‘spect’. A table for the filter response might look like
This is not a particularly realistic table. A proper table should have enough points to ensure that there are a reasonable number of values tabulated over the region of the spectrum to be corrected. With this table in ‘filter.tab’, a filter calibration spectrum can be produced and applied as follows:
Note that ‘idiv’ is used to apply the corrections, which is simply a case of dividing the spectrum to be corrected by the interpolated filter response. The same calibration spectrum may be used for any other spectra that cover the same wavelength range.
A ‘spiketrum’ is a half-way house between a table of values and a spectrum. It has an element for each wavelength value in its range, but only a few of these elements—those that correspond to the table entries—actually have non-zero values. Obviously a spiketrum may be generated very simply given just a table of wavelengths and values at those wavelengths, and a range of wavelength values to be covered. Usually an existing spectrum is used as a template to indicate the wavelength values, and the resulting spiketrum has elements that match those of the template spectrum exactly in wavelength. If a spiketrum is plotted, the result is a set of spikes of varying height—hence the name.
A spiketrum may be turned into a spectrum by interpolating between the spike values.
The ‘gspike’ command will generate a spiketrum from a table of wavelengths and data values, and the ‘interp’ command will interpolate between the points to generate a spectrum. ‘gspike’ also records the nearest values that are just outside the wavelength range covered by the template spectrum, so that ‘interp’ may make use of these as well as the actual spiketrum values. ‘gspike’ also allows some control to be exercised over the data structure of the spiketrum—‘SET’ records included in the table file can cause the UNITS or LABEL objects in the structure to be set to appropriate values. As an example, see any of the supplied flux standard table files (files intended for ‘gspike’ usually have a ‘.tab’ extension); these generally set the units of the spiketrum to match those used in the table, so that ‘gspike’ can produce spiketra whose units are AB magnitudes, micro-Janskys, or whatever, simply depending on what is in the table file used.
For more details about ‘gspike’ tables, see Section 4.5.2.
The main reason for the use of spiketra is that they enable what is essentially tabulated data—instrumental response values at certain wavelengths, as calculated by ‘cspike’, for example—to be manipulated using the standard Figaro routines designed to manipulate spectra.
Since spiketra are really just spectra, they can be plotted using ‘splot’. They may be modified, if necessary, using the fudge commands such as ‘tippex’ or ‘setobj’. If ‘spike’ is a spiketrum and ‘smooth’ is the interpolated spectrum generated from it by ‘interp’, the following sequence will generate a plot of the two superimposed—
Alternative commands to ‘interp’ for interpolating between the points of a spiketrum are ‘spifit’ (which fits a global polynomial to the points) and ‘linterp’ (which uses linear interpolation). The ‘spied’ command is designed to help modify spiketrum points in order to influence the interpolation result—i.e. to fudge the resulting spectrum.
First, a brief word about units. Although it may not always be obvious, the underlying philosophy of Figaro is to provide tools which you may use as you wish, and not to force you to reduce your data in a way imposed by the author of the software. It would be in keeping with this for there to be no restrictions on the units that can be used for flux calibration. However, there are practical limitations.
Magnitude units, such as the AB79 system of Oke & Gunn,3 give results that have a comfortable feel for optical astronomers, and also fall into a numerically convenient range, generally being numbers between 0 and 25. Unfortunately, being logarithmic (not to mention going backwards) they are, while not actually difficult to handle, sufficiently different to linear units to make it impossible to write software that can deal with them other than as a special case.
Rather than do that, Figaro compromises. The flux calibration routines insist on the calibration being performed using linear units. When the flux calibrated spectrum has been produced in such units, it may then be converted to the AB79 scale. This means that rather than have a lot of software that operates on both linear and logarithmic data, we have one routine (‘abconv’) that will perform the conversion.
With that restriction, Figaro will allow you to use whatever linear units you prefer. As will become clearer soon, the units used are determined entirely by entries in the original tables used by the routine ‘gspike’, and if you insist on ‘Joules per minute per square foot per Angstrom’ you are quite at liberty to prepare your own table giving the flux density for your standard at the appropriate wavelengths in these units. You will even find that ‘splot’ will label your axes properly when you plot the data!
As far as linear units go, there are objections to using ‘ergs per sec per square cm per Hz’ and even more so to ‘erg per sec per square cm per Angstrom’ on the grounds that the numbers involved are ridiculously small (and in the latter case can easily go outside the floating point range, which is a serious problem!). So, mainly to be able to avoid having to type commands such as ‘splot high=2.5e-27 low=1.5e-27’, the preferred units for Figaro files are Janskys, mJy, or micro-Janskys.
For most of the flux standards for which Figaro supplies tables, two tables are provided: one in AB magnitudes, usually a direct copy of the published data, and one in (milli- or micro-) Janskys, usually the result of a semi-automatic conversion from the former. If you particularly like erg/s/cm**2/Angstrom, then you can either provide your own flux tables in these units and work from them directly, or you can use the ‘flconv’ command to convert from a spectrum calibrated in Janskys into these units.
The main Figaro data directory (corresponding to environment variable FIGARO_PROG_S) contains a number of files giving the published flux densities of standard stars, all with a ‘.tab’ extension (see Appendix E.1). In addition Jeremy Walsh has made available the Oke and HST standards and copies can be retrieved by anonymous ftp. See Appendix E.2 for details. For example, the file g158m100.tab begins as follows:
The lines beginning with asterisks are treated as comments, and the lines that begin with ‘SET’ are used to set data objects in the file created from this table. An alternative version of this file is ‘g158m100a.tab’, which contains the lines
The functions of the UNITS and LABEL lines should be fairly obvious. Setting MAGFLAG to 1 in a file indicates to ‘splot’ that the data are in magnitude units and so should be plotted with the flux scale reversed. (Note that most ‘.tab’ files actually used by Figaro in fact use ‘.Z.UNITS’ rather than just ‘UNITS’; these were written for the original version of ‘gspike’ where you were allowed to assume that data units always were held in a file in an item called ‘file.Z.UNITS’. Now that Figaro supports NDF format files, this is no longer the case, and the abstract term UNITS is preferred—however, the original format files still work, since the new ‘gspike’ uses a conversion table to handle these explicitly named items.) Tables based on data from, for example, Oke & Gunn’s 1983 paper will also include the line
to indicate the 40 Angstrom bandwidth used by their data. The Filippenko & Greenstein data represent fitted continuum fluxes, so do not have a bandwidth—a point we shall have to return to very shortly. From the ‘.tab’ files already supplied, it should be possible for you to deduce how to create your own, should that be necessary.
There is a Figaro convention regarding the naming of the table files. The ‘ls’ command may be used to find which files are available as tables. The command
will list all the table files supplied in the main directory. Note that not all of these are intended for flux calibration; some may be extinction tables, etc. If in doubt, these are all text files, and should have comments at the start describing their function. So, the command
will list the file for you. You should find that most flux files exist in two incarnations, as implied above, one in a Jansky based unit, one in AB magnitudes. The name of the AB magnitude file ends with ‘a’. So, for example, the files ‘l74546.tab’ and ‘l74546a.tab’ both represent the standard star L745-46A, but the former is in Jansky units—and is therefore probably the one you should use—while the latter is in AB magnitudes. The fact that the name of the object itself ends in ‘A’ is an unfortunate complication that may be misleading.
To anticipate for a moment, the final step in the flux calibration process is carried out by the command ‘spflux’. ‘spflux’ multiplies an observed spectrum by a calibration spectrum, to create a flux-calibrated spectrum. Each element of the calibration spectrum contains a value which is effectively the instrumental response of the detector at a given wavelength, the response being in units of ‘Flux density units per count per second per Angstrom’. The ‘Flux density units’ may be any linear units, e.g. mJy. ‘spflux’ assumes that the spectrum to be calibrated is still in counts, and for each element calculates the wavelength range covered by that element, and then combines that with the counts in that element, the value of the calibration spectrum at that element, and the exposure time for the spectrum, to generate a flux density for the central wavelength of the element. The result is, of course, a spectrum in ‘Flux density units’—whatever they happen to be.
‘spflux’ is straightforward enough; the problem, as with similar functions, is to generate the flux calibration spectrum.
Published standards fall into two distinct classes; they are similar, but differ sufficiently that Figaro provides different sets of functions for dealing with them.
Older, brighter, standards such as those in Oke & Gunn4 are published as tables giving wavelengths and the corresponding flux densities calculated by measuring the observed flux over a range of wavelength centred on the tabulated wavelength. (This is the significance of the .Z.BANDWIDTH value shown in Section 4.5.2.) These tabulated values therefore correspond exactly to what the observed spectrum should look like, even in the presence of absorption features.
More recently, Filippenko & Greenstein5 have published tables for fainter stars where they fit a continuum to the observed data and tabulate the value of this continuum at various wavelength values. These values will not therefore represent the actual observed data in regions where there is absorption, and the concept of a ‘bandwidth’ does not apply.
(You may find it useful to review Section 4.4. before reading on.)
Figaro usually deals with tables, which are tricky to manipulate directly, by turning them into spiketra, which can be manipulated like spectra, although they only have non zero values at points that correspond to those tabulated. No matter which type of published standard is involved, the first step is to generate a spiketrum from the table.
This requires a template spectrum to give the wavelength scale. The one to use is the observed standard spectrum, preferably already scrunched to a linear wavelength scale. (You can use unscrunched data, but it is not advised.) For example, suppose ‘standobs’ is such a spectrum of the standard star HD 84937. There is a table of flux densities for this star (one of the Oke & Gunn standards) in the main Figaro directory, in the file ‘hd84937.tab’. The command
will generate a spiketrum called ‘hdspike’ that can be used by the subsequent steps. If you want a look at it,
will show it as a series of vertical spikes, giving the flux density at each of the tabulated points in the wavelength range of the observation.
‘gspike’ will search directories for the table file in the usual Figaro order: first the default directory, then the user’s Figaro directory (environment variable FIGARO_PROG_U), then the local and Figaro directory (FIGARO_PROG_L), and finally the main Figaro directory (FIGARO_PROG_S). So a table file in any of these will be found, should you need to create your own.
The command ‘cspike’ takes this generated spiketrum (hdspike) and combines it with the observation of the standard (standobs) to generate a new spiketrum, whose values are now the instrumental response sampled at the points of the original spiketrum.
generates this new spiketrum and calls it ‘calspike’. The values in ‘calspike’ are calculated for each point by summing the counts in the observed spectrum over the appropriate wavelength range, dividing them by the wavelength range and the exposure time, and dividing the result into the flux density value given in ‘hdspike’. The units of ‘calspike’ are therefore ‘units per (count per Angstrom per second)’.
‘calspike’ can be turned into the calibration spectrum required by ‘spflux’ by interpolation. The most direct way to do this is just to use the command ‘interp’.
will produce a new spectrum, ‘calib’, which will be the required calibration spectrum. However, you may find that you are not happy with ‘calib’. It may not be smooth enough. It may be that if there are regions of absorption in the spectrum, poor alignment of the wavelength ranges will result in some spurious values. It may seem to be cheating, but the most direct way to get round this is to edit the spiketrum ‘calspike’ using the ‘spied’ command.
will display the spiketrum and allow you to delete points, insert new points, and see the results of spline interpolation or global polynomial fitting to the modified points. When you are happy, you can run ‘interp’ on the modified spiketrum to produce a less honest, but more satisfactory, result.
Global polynomial fitting may seem a trifle crude as a way of interpolating between spiketrum points, but it may be that this gives a better result in some circumstances.
would be an alternative to the use of ‘interp’ to generate the calibration spectrum.
The most direct way to deal with data where the published values represent fitted continuum values is to fit a continuum to your observed data, interpolate the tabulated data, and divide the two resulting spectra to get a calibration spectrum.
Fitting a continuum is a task that is not easily automated. The simplest way in Figaro is to display the spectrum and then use ‘cfit’. For example:
will enable you to use the graphics cursor to indicate continuum points on the displayed spectrum and thus generate the continuum spectrum by spline interpolation. A messy job, but not one you have to do often. Remember that the messy parts of the flux calibration are connected with generating the calibration spectrum, and you only have to do that once; applying it is simple.
The published points represent a smooth curve, so just interpolating directly between them should be quite satisfactory. That is, there should be no need to edit the spiketrum generated by ‘gspike’ prior to using ‘interp’, although ‘spied’ is always available if necessary. This time, since HD 84937 is an Oke & Gunn standard, let’s use G158-100 instead. So the first two steps will be
generates ‘gmfit’ as the interpolated spectrum from the published points.
Dividing ‘gmfit’ by ‘standfit’ will now generate the required calibration spectrum. Note that ‘idiv’ will not do, since the division has to allow for the wavelength range of each element, and for the exposure time. The command ‘caldiv’ must be used instead.
will generate ‘calib’ as the required calibration spectrum.
An alternative way of dealing with this type of data would actually be to treat it as if it were Oke & Gunn data and use ‘cspike’ and ‘interp’ as described in Section 4.5.6). This is probably satisfactory so long as the data does not cover a range with significant absorption features. (If it does, you can always remove the bad points with ‘spied’ before interpolating.)
If you do this, you have to supply a bandwidth, since the table file will not have specified one. You can put one into the spiketrum generated by ‘gspike’ (Section 4.5.5).
will set the bandwidth to 40 Angstroms. The effect of this is going to be similar to smoothing the observed spectrum with a filter 40 Angstroms wide in order to get the continuum values at the points specified by the table. After that, you just go through the ‘cspike’, (‘spied’), ‘interp’ sequence as in Section 4.5.6, finally ending up with a calibration spectrum called ‘calib’.
Given ‘calib’ as the calibration spectrum, no matter how it was generated, it is applied as follows:
which calibrates a spectrum called ‘obs’, creating a resulting spectrum called ‘calobs’.
At the risk of being obvious, you should be aware of what you have created here. Each element of ‘calobs’ is a sample at a particular wavelength of the continuous flux density function; it is not a measure of the total flux within a wavelength bin of finite width, which the original spectra were. Apart from anything else, adding two such spectra has the magic (and totally unreasonable) effect of doubling the magnitude of the object! Generally, this is not something to worry too much about, but it is disconcerting to have a spectrum like this generated with a non-linear wavelength scale. (One has a gut feeling that if the bin covers half the wavelength it should have half the value, and that is not true for these spectra.) That is the reason for the warning about not using unscrunched data.
A calibrated spectrum generated in Jansky units (and others, eventually) can be converted to AB magnitudes by ‘abconv’. Our ‘calobs’ spectrum may be converted by
(‘abconv’ works out the units of the input spectrum by looking at a .UNITS data object in the input file. In some cases, it may not recognise the units—this can happen with spectra that have come from some other system via a translation program. If that happens, and you know that the units are, say, Janskys, you can set the units by hand using
‘abconv’ will recognise ‘mJy’, or anything that contains ‘Jansky’ and the words ‘milli’ or ‘micro’. The same remarks apply to ‘flconv’.)
An interesting test of the system is to calibrate an object using itself, convert the result into AB magnitudes, and then compare the result with a spiketrum generated from the published AB magnitude tables. For example, remembering that our original spectrum of HD 84937 was ‘standobs’, and given that there is a table called ‘hd84937a.tab’ in the main Figaro directory that has the values in AB magnitudes, we can try
The result should be the spectrum of HD84937 in AB magnitudes, with a series of spikes just touching the spectrum, indicating that (at least at the tabulated wavelengths!) the spectrum has been calibrated to the correct value.
A much tougher test would be to use one standard to calibrate another, and then compare the result with the tabulated values. This does of course require that you really do have spectrophotometric data, with no filter changes, with all of the object in the slit in both cases, and with no clouds in the way.
The flux calibration process requires the exposure time of the observation being calibrated. This datum is stored in a standard location in the NDF file holding the spectrum. If your data originates from one of the usual sources, such as the AAT or UKIRT, then the exposure time will probably already be in the correct place in the data structure. However, if your data came from an unusual source then the exposure time may either be present in an unusual location or absent altogether. In these cases you will need to manually add the exposure time to the NDF. The item that has to be added is:
It should be of type _REAL and expressed in units of seconds. Figaro applications ‘creobj’ and ‘setobj’ can be used to respectively create and set the exposure time. E.g. if your data file was called ‘myobs.sdf’ and you wanted to set an exposure time of 90 seconds then you would type:
If the structures ‘MORE’ and ‘FIGARO’ do not already exist in your data file then you will need to create them before creating the exposure time:
For Oke-Gunn type data—
For Filippenko-Greenstein type data—
A number of Figaro functions are available to manipulate complex data, generally with a view to its being used for some process involving Fourier transforms. While there are packaged Figaro routines, such as ‘scross’, which make use of Fourier transforms internally, the functions covered in this section perform the more elementary operations, and can be put together to form a sequence of operations that duplicates the processing performed by, say, ‘scross’, but enabling a finer control to be exercised over the details of the procedure. The general design of this set of routines is based on those provided as part of the SDRSYS system (Straede, 1985). In these notes the term ‘Fourier transform’ is used rather freely; it should be realised that in all cases it is the discrete Fourier transform that is meant.
In Figaro 3.0 the structure of complex data files was changed slightly. Prior to Figaro 3.0, such files contained the modulus of the data as the main data array and held the real and imaginary parts separately. From Figaro 3.0 onwards, the modulus is no longer held as a distinct item, the main data array is the real part of the data, and the imaginary part is still held separately. In practice this should make little difference to the use of these routines, although it does mean that a complex file created by Figaro 2.4 or earlier will not be handled properly by Figaro 3.0. In Figaro 5.1 the format for complex data is changed again slightly. The reason all these changes should make little difference is that complex files are normally intermediate files and are not retained over a long period.
Figaro defines a ‘complex data structure’ in a fairly precise way. A complex data structure is a structure containing two arrays, one holding the real part of the complex data, and one holding the imaginary part. These arrays have to be the same size and shape, although they may have any number of dimensions. The dimensions of the arrays must be such that the Fast Fourier Transform algorithm (Cooley & Tukey, 1965) may be applied to them. Different implementations of the FFT have different restrictions—many require that the number of elements be a power of 2. The NAG library used by Figaro before version 5.0 required that, in each dimension, the number of elements had to have fewer than 20 prime factors, none greater than 19. The routines used from version 5.0 onwards have no restrictions. They are in the Starlink PDA library, the algorithms come originally from the FFTPACK library.
In a Figaro complex data structure the main data array is the real part of the complex data. This means that most of the Figaro routines may be applied to a complex data structure, and if this is done it will be the real part of the data that they operate on. For example, an easy way to look at complex data is with ‘splot’, which will happily plot the real part of the data, quite ignorant of the fact that the data structure in question is complex. To plot the imaginary part, or the modulus of the data, these will first have to be extracted using ‘cmplx2m’ or ‘cmplx2i’. It is important to note that if an ordinary Figaro function is used to change this data, for example by
which multiplies the real array by 2, the imaginary part will be unchanged. Generally, to avoid unnecessary conversions before using the FFT routines, the real and imaginary arrays are held in double precision. However, this is not a strict requirement.
In most cases, one starts with an ‘ordinary’ (i.e. non-complex) file and wants to produce a complex structure in which this is the real part of the complex data and the imaginary part is zero. This is what the command ‘r2cmplx’ does. For example,
will generate a complex structure with the real part taken from ‘myspect’, and the imaginary part will be set to zero.
In some cases, you may want to set the imaginary part of a complex structure as well as the real part. In this case, you can use ‘i2cmplx’, which takes the data from a non-complex structure and uses it to set the imaginary part of an existing complex structure. This means that the complex structure has to be created initially by ‘r2cmplx’, and then the imaginary part can be set by ‘i2cmplx’. So, if you have two spectra called ‘rspect’ and ‘ispect’ which you want to be the real and imaginary parts of a complex spectrum, the sequence is
and the order of operations is important; doing the ‘i2cmplx’ step first will produce a quite different result: the ‘i2cmplx’ will fail, unless it so happens that ‘cmplxspect’ already exists, and even if it succeeds the ‘r2cmplx’ will just produce a new version with a zero imaginary array.
There is no procedure supplied to create a complex structure with the imaginary data taken from a specified file, say ‘ispect’, but with the real part set to zero. It was thought that this would be an unusual requirement. If needed, the following sequence may be used:
The simplest way to ensure that your original data goes smoothly to zero at the ends is to multiply it by a filter that is unity for most of the spectral range but goes down to zero at each end in a smooth manner. The most common form for such a filter is the ‘cosine bell’, and this is what is generated by the Figaro ‘cosbell’ function. (For a detailed discussion, see Brault & White, 1971, and the references they quote).
The only parameter needed by ‘cosbell’ is the percentage of the data that is to be covered by the bell shapes at each end of the data. 10% is a common value to use. ‘cosbell’ uses an input data structure as a template and generates a structure that is the same as the template except for the data itself. Usually, you use the data to which you intend to apply the filter as the template. So, for example, to apply a 10% cosine bell to the data in ‘myspect’,
At present, ‘cosbell’ cannot handle data with more than two dimensions.
Actually taking the Fourier transform of a complex data structure is quite straightforward. The forward transform is performed by the Figaro function ‘fft’ and the reverse transform is performed by ‘bfft’. The only parameters for ‘fft’ and ‘bfft’ are the names of the input and output structures.
For example, to calculate the Fourier transform of the spectrum held in the non-complex structure ‘myspect’, it should be enough to do:
which results in ‘cmplxspect’ containing the ‘fft’ of ‘myspect’. If the power spectrum of ‘myspect’ is required it can be obtained by ‘cmplx2m’. For example, it may be plotted using:
Actually, the modulus is the square root of the power spectrum; the statement that the power spectrum is available represents a slightly cavalier attitude towards the proper usage of terms. However, if the power spectrum is only needed in order to get a feel for the frequency distribution of the original data, then the modulus will generally do just as well.
Often, a power spectrum needs to be plotted logarithmically to produce a sensible plot, so it is quite acceptable to try
The FFT routines generate data with the low frequency terms in the lowest numbered elements in each dimension of the data. ‘fft’ (the Figaro program) re-orders the transformed data so that the zero frequency term is in the centre of each dimension, and the resulting data now goes from -N to +N (where N is the Nyquist frequency). The plots of the modulus made as described above will show a plot with an axis labeled from -1 to +1, (the unit being the Nyquist frequency), that is symmetrical about the centre and which peaks (for most input data) in the centre. New axis structures are created by ‘fft’ to reflect the new axis values. The reason for this re-ordering of the data is that it seems to be easier to visualise complex filters, particularly in 2-dimensions, when the low frequency terms are collected in the middle of the data rather than scattered into the corners.
The reverse FFT is performed by ‘bfft’. Note that the precise definition of ‘forward transform’ and ‘reverse transform’ differ between FFT implementations. The PDA routines, for example, do not introduce any scaling factors between data that is first forward and then reverse transformed and the original data. This means that the sequence
should generate a ‘newspect’ that is exactly the same as the original ‘myspect’, except that any axis information will have been lost.
The previous section sneaked in a reference to the function ‘cmplx2r’, in the hope that what it did was obvious. As the name is intended to imply, ‘cmplx2r’ is the reverse of ‘r2cmplx’, and generates a non-complex data structure whose data array is taken from the real part of the data in a specified complex data structure.
‘cmplx2r’ does not make any changes to the dimensions of the data array.
Analogous to ‘cmplx2r’ are ‘cmplx2i’ and ‘cmplx2m’, which extract the imaginary part of the data and the modulus of the data. Usually, there will be little point in explicitly extracting the real part, since it is already available to most Figaro functions as described earlier. However, the file generated by ‘cmplx2r’ will be smaller than the original file, since it will no longer contain the imaginary array, and this may be a point in its favour.
There isn’t a great deal of point in just transforming data back and forth, unless it’s to get a feel for the errors introduced by such a process. Usually one transforms into the Fourier domain in order to do something useful there. The obvious Fourier domain operations are multiplication (equivalent to a convolution operation in normal space), and multiplication by the complex conjugate (equivalent to a cross correlation in normal space).
Figaro provides a number of functions that operate on complex data to provide a complex result. The operations performed by these should be fairly obvious from their names: ‘cmplxmult’, ‘cmplxdiv’, ‘cmplxadd’, ‘cmplxsub’, and ‘cmplxconj’. The most useful will probably be ‘cmplxconj’, which produces the complex conjugate of a complex data structure, and ‘cmplxmult’, which performs the complex multiplication of two complex data structures.
For example, the cross-correlation function of two spectra may be determined crudely by transforming each, taking the complex conjugate of one, multiplying the two and transforming back. This is crude, because it omits any filtering and preparation of the data, but it will serve as a demonstration of the complex operations involved. If the two spectra in question are ‘spect1’ and ‘spect2’, then
will produce a cross correlation function in ‘corrln’.
Because the data in the Fourier domain is in frequency order, it is often simpler to create a filter in the Fourier domain than to produce the corresponding convolution function in normal space and then transform it. In fact, determining the optimum filter for data is often best done by examining the power spectrum of the data to be filtered and then designing the filter around it. Brault and White (1971) spend some time on this topic.
At present, Figaro provides only one function that generates complex filters, and this is ‘cmplxfilt’. This produces filters of the type described by Hunstead (1980) for use in cross-correlation measurements. These are mid-pass filters formed by taking a Gaussian that peaks at zero frequency and drops to a half height at a specified high cut value, and subtracting a narrower Gaussian that rises from zero reaching its half height at a specified low cut value. That is, the functional form of the data is given by
where u and v determine the low frequency and high frequency cutoff values. From this definition it is clear that u and v are just the sigma values for the two Gaussians, and that what has been rather loosely termed a ‘half height’ is in fact the point at which the Gaussian has a value of exp(-1/2) ˜ 0.6. They are specified for ‘cmplxfilt’ in terms of the Nyquist frequency, so a filter might have, say, u = 0.1 and v = 0.3. The best way to get a feel for these filters is to generate them and plot them superimposed on the power spectrum of the data to be filtered.
If the low cutoff value is specified as zero, ‘cmplxfilt’ generates a low pass filter whose functional form is just
i.e. a Gaussian that drops from unity at zero frequency, having a half width of v. Note that the cyclic nature of the Fourier domain data means that the filter generated is actually symmetrical about the mid point of the data—this is something to be remembered if you want to try to produce other, more specific, filters by generating them yourself. The imaginary parts of the filters generated by ‘cmplxfilt’ are always zero. This means that ‘cmplxmult’ will have the effect of just multiplying both the real and imaginary parts of the data to be filtered by the real part of the filter (obviously, since [a + i b] * [c + i d] = a c + i b c, if d = 0).
‘cmplxfilt’ requires a template complex data structure, which it will use as the basis for the filter it produces. Normally, this will be the data to be filtered, although any data of the same dimensions will do. If the template data is n-dimensional, so will be the resulting filter. The low and high cutoff frequencies will be the same in all dimensions. Since they are specified in terms of the Nyquist frequency, this is usually fairly satisfactory.
So, to take an easy example, an elaborate way of smoothing a spectrum, ‘myspect’ say, by Gaussian convolution, would be
where the 0,0.3 parameters for ‘cmplxfilt’ indicate that a low pass filter is to be created. The 0.3 indicates a cut at around a third of the Nyquist frequency for the data, and the actual value is one best determined by comparison with the power spectrum of ‘myspect’, obtained by performing a ‘splot’ of ‘cmplxspect’ just after the FFT is obtained. Being able to do this is really the main justification for indulging in such a complex procedure when ‘ixsmooth’ could do the same job far faster.
Note that the sequence given above is a pretty complete one. The use of ‘istat’ at the start is to determine the mean level of the data in the spectrum and the result of the ‘icsub’ is to reduce the spectrum to a zero mean level. This reduces the constant component of the power spectrum and should produce a better result. (Note that ‘istat’ sets its output parameter ‘stat_mean’ to the mean value it determines.) The data is apodised by the application of a 10% cosine bell prior to the transform. Both the effects of the cosine bell and the subtraction of the mean level are reversed at the end of the sequence.
As an additional example, consider the following sequence, which attempts to duplicate the functions performed by ‘scross’. It is not an exact duplication, since the internals of ‘scross’ are rather different to those of these FFT routines, and the filters used by ‘scross’ also differ from those generated by ‘cmplxfilt’.
‘result’ now contains the cross-correlation of the two spectra, ‘spect1’ and ‘spect2’. ‘peak’ is a utility produced especially for this sort of sequence, listing the position of the highest peak in the spectrum it is given in terms of a shift from the centre of the first element. This is in fact the relative shift of the two spectra ‘spect1’ and ‘spect2’.
Note that the cross-correlation peak will in fact not be central in ‘result’, but will be at one end or the other. It is often easier to handle data like this if the data is rotated so that the peak is roughly central, and this can be done by the function ‘rotx’.
where 765 is half the number of pixels in the data (here assumed to be 1530). This has the effect of turning ‘result’ into the cross-correlation spectrum that might have been produced by ‘scross’.
Some of the other strange numbers in the sequence may need explanation. The ‘4’s in the ‘sfit’ indicate that the sequence fits a 4th order polynomial to the spectra, and then subtracts away the polynomial fit. This is in fact what ‘scross’ does, unless the fitting is disabled. The 0.1,0.3 parameters in ‘cmplxfilt’ are quite arbitrary, and the optimum values to use are in fact best determined either by looking at the power spectrum of the data to be filtered, or by consideration of the types of features one wants to emphasise in the determination of the correlation (Hunstead, 1980).
The Figaro routines described in this section originated in the need to straighten the distorted spectra produced by image tube detectors. However, even with CCD detectors—which do not suffer the geometrical distortion of the images tubes—some instruments, particularly échelle spectrographs, still produce curved spectra, and the techniques described here can be used to correct these as well. Indeed, nowadays the main application of these routines is to échelle data.
Instruments such as the 2D-Frutti and IPCS that use image intensifiers suffer from various distortions, in particular S-distortion. The effect—and the reason for the name—can be seen clearly by displaying any 2D-Frutti image of a point object, such as a star, on the display. Instead of being a perfectly horizontal line, the spectrum snakes across the image in the shape of a horizontal letter S.
Note that this distortion is actually a two-dimensional distortion—a picture of cartwheel taken through an image tube will show all the spokes bent into S shapes in a radially symmetric manner. However, the difference in pixel scale in the two dimensions for spectral detectors means that—to a first approximation—the distortion can be treated as though it were simply a vertical displacement in the data whose magnitude varies along the spectrum (and to a lesser extent, with position along the slit). In any case, the two dimensional distortion can be corrected by two orthogonal one dimensional corrections: the S-distortion correction described here is one, and the other can be performed as a side effect of a two-dimensional re-binning to a linear wavelength scale. It may be that such full two-dimensional ‘scrunching’ is regarded as overkill—see Section 4.8 for more details.
The process described here is a one-dimensional correction, in which data is re-binned in the direction perpendicular to the dispersion, in such a way as to straighten the spectra.
The tricky part is to determine the exact shape of the distortion, and in Figaro this is done using the command ‘sdist’. ‘sdist’ requires an image with one or more spectra in it. The usual thing to do is to take an observation of a flat-field through a multi-hole mask, in which case the resulting image will show one spectrum for each hole in the mask, each distorted into an S. The spectra will be roughly parallel, but the two-dimensional nature of the distortion will mean that they are not exactly so; this is the reason for using multiple spectra: it allows the change in distortion with slit position to be determined. The ‘cdist’ command can then be used to apply the results of the distortion analysis performed by ‘sdist’ in order to correct an image.
The usual sequence of operations is as follows: assume ‘holes’ is a multi-hole calibration image, and ‘object’ is an image to be corrected.
There are a couple of decisions that have to be made in this sequence. Firstly, there is the question of just how many of the hole spectra should be used. Generally the answer is that all should be used, unless they are going to cause problems—a spectrum that gets too close to the edge of the image will probably not be traced properly, for example. Then there is the use of the ‘columns’ parameter in ‘sdist’ to help control the tracing of the data. In most cases, the main point about the ‘columns’ parameter is that making it larger reduces the number of points to which the polynomials are fitted, and so speeds up the process considerably. However, if the spectra are particularly weak, increasing ‘columns’ will also improve the signal to noise of the cross-section profile of the spectra. (‘sdist’ sums the data over the specified number of columns, then tries to locate the centre of the resulting profile, assuming that the best starting guess for the centre is the value obtained from the previous fit.)
‘sdist’ can use a number of different algorithms for tracing the spectra, selected by the ‘trace’ parameter. The original algorithm used by ‘sdist’ is the G(aussian) option, and this is most suitable if the profile of the summed data is roughly Gaussian (as it usually is for stars observed using a slit). If the profile is roughly ‘top hat’ in profile (if a dekker has been used, for example), then either E(dge) or C(enter of gravity) will probably be better. Edge fitting tends to produce rather jumpy fits, particularly if the edges are very sharp and therefore cannot be located to better than a pixel. The centre of gravity fit (which takes the COG of the data within the edges) is usually smoother, but can be inaccurate if the top of the data is not flat. A ‘ystract’ through the data, followed by a ‘splot’ can give a good feel for the data profiles, and the diagnostic display produced by selecting the ‘softd’ option can be very helpful here. The information produced by specifying the ‘diagnostic’ keyword is really for debugging new trace algorithms and is unlikely to be of general use.
In the end, there is no substitute for watching the results on the display. Remember that the polynomial will be unconstrained outside the range where it could trace the spectrum, and this may well show up in the final results from ‘cdist’. So see if the spectra seem to have been fitted acceptably, and if not, either ignore the bad ones or try to fine tune the fit using the ‘sdist’ parameters. The display that ‘sdist’ can produce on a line graphics device is a useful diagnostic. (In Figaro 3.0 ‘sdist’ could overlay the fit on a previous display done with ‘image’. Sadly this option is not available in the current version.)
‘offdist’ may be useful if there is a slight linear offset in Y between the calibration data and the data to be corrected. This is not usually important if the data is to be corrected by ‘cdist’, but other routines such as ‘maskext’ can make use of ‘sdist’ results in a way that makes such offsets (often the result of guiding changes) important. ‘offdist’ adds an offset in Y to the results produced by ‘sdist’.
The final step, the correction of the object data, is an automatic process. So if a number of images are to be corrected there is a lot to be said for doing this as a batch job.
Image tubes are not always tremendously stable; it may not be reasonable to assume that the details of the distortion do not change with time. This means that if there is only one multi-hole calibration image, probably taken at the start of the night, a distortion analysis based on this single early image may become progressively less correct as images from later and later in the night are processed. And, of course, there is always the possibility that no such image was taken at all.
In these cases, it is always possible to produce a distortion analysis from a spectrum of a single point object. An image of such an object can be used to calibrate itself, from a distortion point of view. The disadvantage of this over the use of a multi-hole exposure is that the correction based on a single object will be exactly correct only at the position of that object in the slit—the sky nearby will not be quite so well corrected.
So the choice of using objects to calibrate themselves, against the use of a multi-hole calibration image, is a trade-off of one source of error against another. The choice has to be yours.
The operation performed by ‘cdist’ involves working on the image column by column. On a virtual memory machine, this is to commit the cardinal sin of accessing an array in the wrong order—against the grain of the virtual memory system. When large images are processed, the result will be excessive page faulting by the process, together with a dramatic increase in the time taken (both CPU time and elapsed time go up), and an even more dramatic increase in the anger levels of the other users of the machine, since the excessive paging will begin to saturate the disk I/O system. This condition is described as ‘thrashing’.
To reduce this faulting, it is possible to rotate the image before it is processed and rotate it back afterwards. ‘cdist’ has a hidden parameter ‘rotate’ which makes it do this automatically. The final results obtained are the same whether or not ‘rotate’ is specified, but the efficiency of the operation can be quite different. Specifying ‘rotate’ adds the overheads of the two rotations, but makes the correction work properly with the virtual memory system. For small images, where the correction algorithm probably does not induce thrashing anyway, use of pre and post rotation will be inefficient since it adds the overheads of the two rotations. However, these overheads are small for small images, and it is probably acceptable to always specify ‘rotate’. This is why ‘rotate’ is a hidden parameter, true by default.
Just what constitutes a ‘small’ or a ‘large’ image—i.e. the size of image at which one starts to gain by specifying ‘rotate’—is hard to say. The best guide is trial and error. As it usually is.
Firstly, if you cannot display grey-scale or colour, you still can use ‘icont’ on any line graphics terminal. Its display is equivalent to that of ‘igrey’ and can be used by ‘igcur’ (not ‘icur’).
It is—just—possible to make do without any display at all, should one not be available. ‘sdist’ picks up the number of spectra and the pixel coordinates from global parameters ‘npixels’, ‘xpixels’ and ‘ypixels’. You can create these with ‘creobj’, or delete them with ‘delobj’ if necessary. ‘xpixels’ and ‘ypixels’ are vectors, if these are too short, delete and re-create them. Once created you can assign values to these parameters using ‘setobj’:
(If ‘ADAM_USER’ is not set, use ‘$HOME/adam’ instead.)
Simple statements such as ‘the wavelength of channel n is so many Angstroms’ are not as unambiguous as they appear. This statement is really a slightly simplified version of ‘the wavelength range covered by channel n is from so many Angstroms to so many Angstroms’ but some precision has been lost in the simplification. It presumably means ‘the wavelength of the centre of the range covered by channel n is so many Angstroms’—or does it mean that ‘the wavelength of the start of the range covered by channel n is so many Angstroms’?
These notes are intended to explain the conventions used by Figaro. Please note that they are all somewhat arbitrary and in some cases are historical rather than ideal. However, they are consistent and they are precisely defined. Hopefully, these notes will also serve as a brief guide to the use of Figaro for wavelength calibrations. If they seem a trifle pedantic, this is because experience has shown that the whole subject is a real can of worms.
Arc-fitting is also available with Specdre (Section B.7) which is now part of Figaro. It is reputed to be far superior to that described in this section. It is currently only useable with Specdre commands but not with other Figaro commands.
Another alternative for arc-fitting is Molly, a semi-interactive wavelength calibration program written by Tom Marsh. It allows the first of a series of similar arcs to be calibrated interactively, then calibrates the rest automatically. In the right circumstances it can be a great time-saver.
Most of the problems arise because
The sad fact is that there is no obviously correct method of mapping the continuous channel numbers (‘x’) onto the integer channel numbers (‘n’). Mainly it’s a question of selecting the zero point. The Figaro convention is that
This all sounds mind-blowingly trivial, but it is possible to produce other conventions. In particular, if the left edge of channel 1 were taken as the x = 0.0 point, the centre of channel 200 would be at x = 199.5. This is in fact the convention for the NDF format, used by Figaro. It is of little practical consequence: when pixel coordinates matter, they are never implied from pixel numbers but stored explicitly. The n-th stored value always applies to the centre of the n-th pixel.
Figaro tries to hold wavelength information in as general a way as possible. If a file called, say, ‘spect’ contains a wavelength-calibrated spectrum, then there will be two data arrays in the structure held in the file. The main data array will be an array containing the actual spectrum, and the X-axis array will be an array of related size containing the wavelength information.
Each element of the X-axis array contains the wavelength value for the corresponding element of the spectrum—that is, for the centre of the corresponding data bin. For uncalibrated data, the X-axis array will usually not exist.
The ‘hdstrace’ command may be used to examine elements of either the data or wavelength arrays, and the ‘ilist’ command can be used to tabulate them.
A single arc spectrum can be generated by summing successive rows of an arc image, using the ‘extract’ command. The wavelength to channel number relationship can then be determined using the ‘arc’ command. This invokes an interactive line identification program which gets the user to select lines in the arc and to specify their wavelengths. ‘arc’ makes full use of interactive graphics and a running fit to make the manual identification process as easy as possible, and it also has an automatic line identification feature that can be used to add lines to a fit. ‘arc’ is described in more detail in Section 4.10.
Once a satisfactory fit has been obtained, ‘arc’ can be used to set the wavelength scale in the arc spectrum being identified. This is the importance of the question ‘make use of this fit?’ that the user is asked at the end of ‘arc’. Once the wavelength scale has been set, subsequent plots of that arc spectrum made by ‘splot’ will be made with a scale in wavelength units.
Note that Figaro programs are all designed to work with data whose X scale either increases or decreases. However, as a practical point, it is regarded as more usual to arrange data so that wavelength increases with pixel number, and there may be bugs still lurking in some routines that will only turn up when used with reversed data. You are therefore strongly recommended to use increasing X scales, if necessary using the ‘irevx’ command to reverse spectra and images.
If ‘arc’ is used to fit a spectrum very similar to the one previously fitted, it is enough to make use of the previous fit (by answering ‘yes’ to the ‘previous’ prompt) and then re-identify a line used in the previous fit. ‘arc’ will ask if the previous identification is to be ignored, or the new identification, and as a final option will assume that there is a shift in the data relative to the identifications and will re-analyse the line centres accordingly. This makes the analysis of a sequence of arcs quite simple.
Of course, the point of getting such a fit is to apply it to other spectra. The command ‘xcopy’ copies the wavelength scale from one spectrum to another. Normally the spectrum whose wavelength scale is copied is a fitted arc and the spectrum to which it is copied is that of an observed object.
Typically, an observation of an object is bracketed by arc observations. These different arcs will inevitably give slightly different fits, either just because of random variations, or because of some instrumental drift. Particularly if there is drift involved, it may be better to use ‘xcopi’ rather than ‘xcopy’. ‘xcopi’ takes as input two spectra with wavelength scales and generates a new scale that is an interpolation between them. It then applies that interpolated scale to the object spectrum.
It is possible to rebin data whose wavelength/channel number relation is known in such a way that the relation is linear. Such an operation is known as ‘scrunching’, and is performed by the ‘scrunch’ command. Scrunched data is generally easier to handle than is un-scrunched. For example, if spectra with slightly differing wavelength scales must be added together they should be scrunched to a common wavelength scale first.
Similarly, there are cases where one needs data rebinned on a logarithmic wavelength scale—cross-correlations to determine redshift need this sort of scale, for example—and this can also be performed using ‘scrunch’.
‘scrunch’ has parameters that specify the start and end wavelengths for the resulting scrunched spectrum, and the number of elements (bins) it should have. Note that the wavelengths are those of the ‘centres’ of the end bins; it is easy to miscalculate the number of bins by 1 when aiming to get a linear dispersion of exactly so many Angstroms per bin. The target wavelength range may also be specified in terms of the start wavelength and the increment, by means of the rather crude convention that if the end wavelength specified is less than the starting wavelength, then that ‘end wavelength’ value will be taken as the wavelength increment. (This behaviour may be controlled explicitly by two keywords ‘final’ and ‘increment’, which are used to specify the interpretation to be placed on the ‘end wavelength’ value.)
‘scrunch’ has one keyword that can confuse users. If the spectrum being scrunched is in flux units, then the total flux in the spectrum should be conserved by the scrunching process. Note that the distinction is between flux units and flux density units, so ‘raw counts’ are a flux unit. If a spectrum is in flux density units, then scrunching should not change the absolute value of the data. If the scrunching is simply such that it doubles the number of bins in the spectrum, then an element that has a value of, say, 10 mJy should become two elements each with a value of 10 mJy, since Janskys are a flux density unit. Conversely, if the spectrum were still in raw counts, then an element with 10 counts should become two elements each with 5 counts. The default for the ‘flux’ keyword in ‘scrunch’—whose prompt is the perhaps confusing ‘conserve flux (as opposed to mean counts)?’ is ‘yes’, and this is correct for raw counts. To confuse matters further, there is a concurrent parameter ‘mean’ with the opposite meaning to ‘flux’.
The previous sections have all been concerned with single spectra, and this is usually satisfactory where images can easily be reduced to spectra by use of ‘extract’—i.e. where there is no significant change in wavelength/channel relationship across the subset of the image that is being collapsed into a single spectrum. However, there is an alternative approach, which is to scrunch each cross section of an image separately. This may be needed for data on extended objects, or it may simply be regarded as a better (more strictly correct) means of proceeding.
To do this, one needs a good fit to every cross-section of an arc image. This is performed automatically by ‘iarc’, starting off from a manual fit produced by ‘arc’. The usual technique is to extract a spectrum from the centre of the image, probably adding together a few cross-sections to get good signal to noise, and then fit that using ‘arc’. It is important to have as many lines as possible in the fit; it is even worth adding lines in any sparse areas to help lock down the fit, even if their identifications are uncertain—just treat the fitted values as exact. ‘iarc’ does have an automatic search for such ‘lock’ lines, but it is usually better to pick good ones manually.
By default, ‘iarc’ performs its analysis using the same line width parameter, ‘sigma’, as was used in the original fit performed by ‘arc’. However, the algorithm used by ‘iarc’ is quite sensitive to the value of this parameter, and better fits may sometimes be obtained by varying it. For this reason, ‘iarc’ has a hidden parameter ‘rsigma’ which may be used to over-ride the ‘arc’ value. If you have problems with ‘iarc’, there is a hidden keyword, ‘detail’, which causes the details of the fit for each cross-section to be put out. This can be used to see which lines the program is having trouble with. ‘iarc’ will drop a line from its search list if it fails to find it in a a given number of successive cross-sections—this is to prevent it accidentally picking up the wrong line due to distortion moving another line towards the position where it lost the original line. The number of cross-sections is given by the hidden parameter ‘gap’—normally, ‘gap’ is set to 1, meaning that the second time in succession a line is missed, it will be dropped. Setting ‘gap’ sufficiently high will effectively disable this feature—this may be the best way to deal with data where there is little distortion, but where a number of cross-sections may have poor data; fibre data can be of this type. If ‘gap’ is set to zero, a line will be dropped immediately if it cannot be found in a spectrum. In order to pick up as many lines as possible in each cross-section, ‘iarc’ uses a technique known as ‘spreading’—it first looks for the lines using a larger sigma value than that specified, and then refines the fits using the specified sigma value. In some cases this is undesirable, and may be bypassed by the use of the hidden keyword ‘nospread’.
Sets of spectra which have been taken in some long slit mode with a distorting detector, such as an image tube, will differ from spectrum to spectrum but will do so smoothly. Sets of spectra taken with a fibre system may vary discontinuously, with sudden linear shifts between spectra. The ‘xcorr’ keyword in ‘iarc’ is designed to be used with such discontinuous spectra. It causes ‘iarc’ to attempt to determine a linear shift between successive spectra by cross-correlation and to apply that shift to its expected arc line positions before searching for them in the new spectrum. ‘xcorr’ should be used for fibre and similar data, but is probably inappropriate for other data.
Given a fit from ‘iarc’, ‘iscrunch’ can be used to scrunch an image of an object. It is actually a good idea to scrunch the original arc and see what the results look like—lines waving about in some regions show that the fit is badly constrained in those regions. ‘iarc’ only performs a sequence of individual one-dimensional arc fits, rather than a two-dimensional fit to the whole image at once, so there is little constraint that the fits be continuous from one cross section to the next, other than that imposed by the lines themselves.
Note that the ‘iarc’ results can only be used by ‘iscrunch’; they are not written into the data structure in a way that can be used by other routines such as ‘splot’.
For object images bracketed by arc images, ‘iscruni’ can be used instead of ‘iscrunch’. This uses wavelength values obtained by interpolation between two different ‘iarc’ fits, in an similar way to ‘xcopi’ as described earlier.
For individual spectra:
For images:
The ‘extin’ command is used to correct a spectrum for atmospheric extinction. It requires a ‘coefficient’ spectrum—a spectrum giving the atmospheric extinction coefficient for each element of the spectrum. This coefficient spectrum is normally generated by interpolating a spiketrum formed from a suitable table of extinction coefficients.
One such table is ‘extin’, (actually the file ‘extin.tab’ in the FIGARO_PROG_S directory). This contains the coefficients given for Palomar by Hayes & Latham (1975). So for Palomar data, the following sequence will perform the extinction correction:
generating a corrected spectrum, ‘cobject’, from the original uncorrected spectrum, ‘object’. The first two steps do not need to be repeated for subsequent spectra, so long as the wavelength range covered remains unchanged.
(An alternative extinction table is ‘palomar’, which contains the coefficients used by TYB’s Forth system. This differs slightly from ‘extin’, particularly around 9000 Angstroms, where it attempts to correct for some atmospheric features ignored by Hayes & Latham. The differences can be seen by generating and then plotting the spiketra produced by the two tables.)
The ‘aaoext’ table gives the standard AAO extinction coefficients—those used by the SDRSYS data reduction system.
Hayes, D.S., Latham, D.W.: 1975, Ap.J. 197, 593
This section describes the use of the ‘arc’ program in Figaro. ‘arc’ is an interactive arc fitter, which displays an arc spectrum on the current soft graphics device and gets the user to indicate arc lines and enter their wavelengths. Given sufficient lines, a polynomial fit of wavelength against channel number can be performed, and the results of this fit used to fill the wavelength array (the .AXIS(1).DATA_ARRAY). Section 4.8 describes the way ‘arc’ fits can be used.
‘arc’ is primarily designed as an interactive arc fitter, but it does have an an automatic line finding capability. This description will emphasise the interactive aspects of ‘arc’, however, since this automatic capability is intended to help add lines to an already good fit, and so the first and most important thing is to get a good fit manually.
‘arc’ has a number of features that are intended to make it particularly simple to use. These are not always obvious to the user who simply types the command ‘arc’ and waits to see what happens next; hence this description.
To get the best out of ‘arc’ you need a comprehensive line list for the arc you are trying to identify. Arc line lists are text files containing lists of arc line wavelengths. They have the extension ‘.arc’ and are usually held in the main Figaro directory (FIGARO_PROG_S; see Section 3.4). There may also be files in the local Figaro directory (FIGARO_PROG_L), or you, the user, may have your own in your user Figaro directory (FIGARO_PROG_U), or in the default directory. As an example, the file ‘argon.arc’ in the main Figaro directory, begins:
and carries on in the same vein for quite some time. The format is very simple—each line of the file has a wavelength, in free-format, and blank lines and lines beginning ‘*’ are ignored. Note that no line strength information is used. The importance of this line list will become clearer as we go on.
When the ‘arc’ command is given, one of the first prompts asks ‘Type of arc?’. ‘arc’ wants to know which arc line lists to read in. You can specify up to three arc types, although it is unusual to use more than one. If your reply is ’arc1,arc2,arc3’, then ‘arc’ will look in the various Figaro directories for the files ‘arc1.arc’, ‘arc2.arc’ and ‘arc3.arc’. It will then read in the arc lines from all three. The idea here is that if you have, say, a copper-argon arc, you might give the type as ‘copper,argon’, and ‘arc’ can make use of two separate files, one for each element. In practice, it would be better to have a single file called ‘cuar.arc’, set up specially for the arc you are using, and to reply ‘cuar’. If you have no suitable arc file, reply ’NONE’.
‘arc’ will also prompt for the arc line half width in pixels (the parameter called ‘sigma’), which it uses as a guide when finding the centres of the arc lines, and for the initial order of polynomial to use for the fit. ‘arc’ performs running fits, and unless you are using a previous line list (see below), the first fit will have to be made to a low order, simply because there will not be enough lines for a higher order. Once enough lines have been identified, ‘arc’ will start to use the order you specified initially. You will be able to change the values for both ‘sigma’ and ‘order’ interactively during the fitting process.
‘arc’ will also ask if lines from the previous fit are to be used. During a fit, ‘arc’ keeps writing out to disk the positions and wavelengths of the lines identified so far. (‘arc’ always writes to the file ‘arlines.lis’, thus any such file should be deleted or renamed before you run ‘arc’.) If you reply ‘yes’ to this prompt, ‘arc’ will read in the file giving the previous list of identified lines. This allows you to start again where you left off the previous time—either because ‘arc’ or the computer crashed (perish the thought!) during the previous fit, or simply because, on reflection, you are unhappy with the fit you obtained and want to experiment with other orders, a different selection, etc. The default for this previous list file is always ‘arlines.lis’, since that is the name of the file that ‘arc’ writes. However, you may specify a different file. You may, for example, have one file that you want to use as the basis for fits to a number of arcs, which you have renamed from ‘arlines.lis’ (the name ‘arc’ will have given it originally when it was produced) to, say, ‘basefit.lis’.
You can also use a previous fit if the arc you are identifying is very similar to the previous arc, but is shifted slightly. This is often the case with arcs taken at intervals throughout a night. ‘arc’ will notice if the file that was used to create the previous fit is different to the file you are fitting. If this is the case, you are prompted for the ‘xcorr’ keyword, which allows you to request that ‘arc’ locate the previous spectrum used and attempt to determine a linear shift by cross-correlating the two arcs and applying the determined shift to the arc line list. (If you are going to use this option, specifying a slightly larger ‘sigma’ than that used for the previous arc will give the line finding algorithm a little more flexibility when it looks for the lines at their shifted positions.) This is a simple operation, but the shift is determined over the whole of the arc. There is an alternative, messier, way to indicate a shift to ‘arc’ on the basis of a single selected line, but this is described later.
‘arc’ will read in the specified spectrum and display a portion of it on the graphics device. Initially, the portion will be 200 channels long; you can change this should you want to. You will be invited to use the cursor to indicate an arc line.
Normally, you move the cursor until it is close to a line whose wavelength you know (you will often find it useful to have a hard plot of the whole arc in front of you as you perform the fit) and select it by hitting the space bar. (Strictly, you can use any key that does not have some specific function, but there are rather a lot that do and the space bar is a safe one to use.) ‘arc’ will then try to find a line close to the point you have indicated, and if it finds one will show you its centre on the display with an arrow. The algorithm used to find a line centre is one described as ‘convolution with the derivative of a Gaussian’, and it incorporates some requirements as to just what constitutes a line, which can lead on rare occasions to its being unable to find a line centred near the position indicated. You can think of the algorithm as a fit to a Gaussian of fixed sigma.
You will be told the channel number of the line, and asked to enter a wavelength. Strictly, you enter a wavelength followed by a type, e.g.
The type should be one of the names specified in response to the ‘What type of Arc’ prompt. ‘arc’ then looks in the table for that type and selects the line whose wavelength is nearest to the one you specified. If you do not specify a type (which is by far the most usual case), ‘arc’ uses the first of the types. Since in most cases only one type was specified, this means that ‘arc’ will normally just look in the one table that it read in. Which is what you would expect. So the response is usually just a number, e.g.
‘arc’ will then tell you what line it assumes you mean:
If you reply ‘y’, ‘yes’, or just hit the return key, ‘arc’ will use that wavelength. If you reply ‘no’, or ‘n’, it will ask you for the wavelength again. If you just hit return in response to the wavelength prompt, the line will be deleted and you will be back with the cursor looking for another line.
There are some important alternative ways of specifying the wavelength of the line.
‘e’ being short for ‘exact’—‘arc’ will use that value as the wavelength. If you have specified ’NONE’ as the arc type, ‘arc’ will assume that all wavelengths given are exact values.
What this means in practice, is that with a good line list, you can identify enough lines to tie down the fit fairly well (it helps to do some lines at one end, then at the other, then work in), and then just respond ‘f’ to each wavelength prompt. ‘arc’ will then use the line in the line list closest to the fitted wavelength. This simplifies the process considerably, but it is still under your control, and you can intervene if a new fit shows a considerable increase in RMS.
When you are asked to indicate a line using the cursor, you can hit—instead of the space bar, which just selects a line—any one of a number of keys that have specific functions.
For example, ‘n’ moves you on to the Next portion of the arc. ‘b’ moves you Back to the previous section. ‘m’ Moves you to a specific X value. (Note that the value is not always in channels—if the spectrum already has wavelength data associated with it, the X values will be in Angstroms. This is sometimes a useful feature, sometimes an irritating one.)
This is how you move around from one section to another. ‘l’ lets you change the Length of the portion displayed at any one time.
‘d’ deletes the identified line nearest the cursor. If you see the RMS for the fit shoot up, hitting ‘d’ will delete the last line identified—so long as you haven’t moved the cursor since you selected it.
An important command is ‘q’, which quits the identification process and moves on to the fitting and line editing stage. Note that this does not commit you to anything—you can always come back to the interactive selection stage.
If the normal line centring algorithm cannot find a line, you can use ‘c’ to force a line identification where the line centre is determined by a centre-of-gravity analysis. This is a time-consuming operation, and is not really recommended for normal use.
Strongly not recommended is the use of ‘e’, which allows you to indicate the line centre yourself on an Expanded plot. This command does nothing clever at all—it just takes the position you indicate.
Like most Figaro programs that involve key-selected options, ‘arc’ will respond to ‘h’ or ‘?’ by printing a list of the options it accepts.
The ‘q’ command takes you to a more conventional stage where the fit is repeated and the results displayed. As soon as you hit ‘q’ and confirm that you indeed want to move to this next stage, a fit is performed to the lines identified so far and the results are displayed. For each line, the calculated and actual wavelengths are given, together with the discrepancy in Angstroms, and the RMS for the fit that would be obtained if that line were deleted from the list of lines. Finally, the RMS for the current fit is displayed. Remember that the RMS is in terms of Angstroms.
If the fit is bad, a look at the table of residuals output may help to indicate which identifications are at fault. The ‘RMS without this line’ is a particularly useful figure when there are relatively few lines in the fit; it saves you the effort of making a tentative deletion and a refit to see if the result is really an improvement. With a large number of lines, the RMS is less dependent on any one line, and this figure becomes less useful.
You are given a number of options at this point, all of which are listed in a single prompt. This is described as the ‘menu’ prompt. The options are as follows:
The first letter of each command is sufficient.
When you select the ‘Auto’ option in the menu, ‘arc’ tries to use the fit you have obtained already as a starting point and tries to find lines in the arc that match lines in the tables. The algorithm used is very simple, and is based on the principle that the automatic fit should not add lines that will make the fit significantly worse than do the lines you have already got—and are presumably happy with. There are only two parameters (at present) involved and, really, only one is important.
‘arc’ takes each pixel in the spectrum in turn. If that pixel is more than ‘chfact’ times the current sigma value from any line already found, it uses that pixel as the starting point for a line search. This is exactly as if you had selected that pixel with the cursor during the interactive part of the process. If anything resembling a line can be found, it calculates its wavelength and looks in the line tables for a line close to that wavelength.
A line is accepted if the discrepancy between its calculated and tabulated wavelength is less than ‘sigfact’ times the current RMS value. It is this that means that the criterion for accepting new lines is based on how their wavelength discrepancies compare with those for the lines that have already been accepted. ‘sigfact’ is the more important parameter of the two. The default value of 3.0 means that the automatic search can make the overall RMS of the fit somewhat worse, but it will give the program a fair go at finding some lines. Setting ‘sigfact’ to 1 or less, which you may do with the ‘Modify’ menu option, ensures that the automatic search will not make the fit worse, but it will probably not find many lines either.
Just how best to use the automatic line finder is a matter of experience and, probably, opinion. At the moment, it is a relatively new feature and so the experience is lacking, even if the opinion is not. It does, however, seem fair to say that the better the original fit, the more likely the automatic fit is to make correct rather than misleading identifications. The ‘Xauto’ menu option, which causes all automatic line identifications to be deleted, at least means that you can experiment a little without doing irreparable damage to your fit.
One approach is to let the automatic fit loose, and then tidy up after it, deleting those lines that it found but that you don’t like the look of. You can do this from the line list produced by the ‘Fit’ menu option, using the ‘Edit’ option to delete the lines affecting the fit the most. Lines found automatically are flagged in the list by a plus sign (and are shown in ‘arlines.lis’ with a (A) symbol).
An alternative, and quite a simple operation, is to return to the interactive selection (the ‘Reselect’ menu option) and examine the lines found. The lines found by the Auto option are displayed with their wavelengths in parentheses, so it is fairly straightforward to run through the spectrum with the cursor, hitting ‘d’ at every bracketed line that looks like nothing more than a noise bump.
If you do not return to the interactive stage, you will be asked if you want to make use of this fit. If you reply ‘yes’, ‘arc’ will generate a set of wavelength values for the spectrum and either create a new output spectrum with these values in the X data array, or simply put these values in the X array of the current arc spectrum. See Section 4.8 for details of how you use these values. If you are going on to ‘iarc’, you do not need to do this. If you are going on to ‘scrunch’, you do need a spectrum with an array of X wavelengths.
Finally, you have the option of producing a hard plot showing the identifications you have made. You can also produce a hard-copy of the dispersion curve—the plot produced by the ‘Disp’ menu option. Note that if you produce both plots, they will be produced as separate files, both of which will have to be output to the hard copy device.
If you identify one arc, and then move on to another arc which is essentially the same except for a linear shift, you can make use of the cross-correlation option (the ‘xcorr’ keyword described earlier), or you can try the following sequence, which gives you a little more control over the way the shift is determined:
This section is more or less the documentation by J.G. Robertson, dated 2 February 1987 as included as on-line document in Figaro 3.0.
‘abline’ is an interactive program which runs within the Figaro data reduction system, and whose main purpose is to find wavelengths and equivalent widths of absorption lines. The facility for fitting a polynomial to the continuum may also be of use in other situations. As well as wavelengths and equivalent widths the program estimates width and asymmetry parameters for each line.
The program does not fit Gaussian profiles to lines; in fact this approach is specifically avoided, since it results in model-dependent wavelength and equivalent width values. In practice one has no reason to think lines are even approximately Gaussian, e.g. the Voigt profiles of saturated or damped lines, or the multiple component lines typical of QSOs.
The wavelength estimator used is the median (i.e. the value that equally splits the area of the line). Not only is this more resistant to the effects of low level wings on the line than is the centroid (mean; centre of gravity) but also it can be shown that the median is actually a statistically better (less noisy) estimator than the centroid for locating absorption features. For further discussion of the estimators used for all the parameters see Publ. Astron. Soc. Pacific 98, 1220, (1986).
The program can also be used on emission lines, and will find median wavelength, equivalent width, line width and asymmetry. This version of the program does not calculate uncertainties of the wavelength or equivalent width.
Before starting ‘abline’, make sure the soft and hard plot devices are appropriately specified, using the Figaro commands ‘soft’ and ‘hard’.
When you are in Figaro, just type ‘abline’ to run the program. Only interactive use (not batch) is feasible. ‘abline’ reads data from one-dimensional spectra only. There are three options for handling the continuum fit in the vicinity of the required line:
‘abline’ needs a graphics terminal for soft plots and cursor input (as specified by Figaro command ‘soft’). It will be useful to have a hard copy plot of the spectrum on hand before entering ‘abline’, so that you can see the approximate wavelengths of lines to be analysed, and reasonable values for the wavelength range over which to fit the continuum. ‘abline’ assumes that the .AXIS(1).DATA_ARRAY structure of the input spectrum gives the wavelength in Angstroms, although other units (e.g. velocity) will work if one suitably re-interprets the wavelength, equivalent width and line width outputs. The .AXIS(1).DATA_ARRAY structure must be linear (i.e. scrunched) because this is assumed in the calculations. Conversion of the .AXIS(1).DATA_ARRAY structure to vacuum heliocentric or LSR can be made using ‘vachel’ (and ‘scrunch’) if necessary, before using ‘abline’. For a noisy spectrum it may be useful to slightly smooth the data used as input to ‘abline’ (this does not bias the wavelength or equivalent width; it obviously does affect the width parameter and to a lesser extent the asymmetry).
‘abline’ has a number of parameters. These may be set in the command line, since they are ordinary Figaro parameters and keywords, but since ‘abline’ is inherently an interactive program it is expected that most users will just type ‘abline’ and allow themselves to be prompted by the program. The prompts in this case are as follows—
If the continuum is to be fitted in this run, the system next asks:
Alternatively, if you elect to read in a pre-computed continuum, the program prompts
At this stage the program examines the .AXIS(1) structure of the input spectrum, and displays its label and units for the user to check (e.g. they will usually be ‘Wavelength’ and ‘Angstroms’). The first and last values of the .AXIS(1).DATA_ARRAY structure are also displayed for checking.
The next parameter requested from the user controls how the wavelength bounds of the absorption line will be found:
The next parameter requested is
You are now asked whether you want to save the continuum fit(s) to a separate spectrum, and if so, what the name of this spectrum should be. Normally one does not want to use the continuum after the present ‘abline’ run; in this case answer ‘no’.
You now reach the loop in which individual lines (all from the same input spectrum) are analysed. A help facility explains the available commands. The program prompts with
The most usual response is a floating point number: This is the approximate wavelength of the next line to be analysed. After a further question to find the name of this line (which will label the hard and soft plots) a soft plot will be displayed. Its wavelength range will be as specified above for ‘wavelength range...’ and it will be centred at the wavelength specified here, unless this is too close to one end of the spectrum, in which case the plot will extend just to the end of the data. Further progress from this point is discussed below.
Other responses enable you to change the way the program fits the continuum or analyses the line. It is useful to be able to change these parameters here because different lines may require different treatment, and we do not want to have to re-enter the program to specify these. The other possible responses are (where xxx is a floating point number and nnn is an integer):
The changes made will remain in force for all subsequent lines analysed unless they are changed again.
At this stage in analysing a line you have a soft plot of the selected region. You now have to use the cursor to specify the exact channel ranges to be used for fitting the continuum (unless a pre-computed continuum fit is being used). It is clearly necessary that the continuum fit does not include the absorption line itself; there may also be other absorption lines or defects nearby which you want to explicitly exclude from influencing the continuum fit. This is handled by allowing you to specify up to 10 wavelength sub-segments (within the plotted range) which will be used as input to the continuum fit. Type ‘cont’ and in answer to the rather extensive prompt ‘Selection of sub-segments for continuum fitting ....’ Do the following:
Having finished sub-segment selection, a polynomial fit is made to the specified continuum channels. Channels whose data values deviate from the fitted curve by more than the ‘multiple of sigma for continuum point rejection’ are rejected and the fit recomputed. This process is repeated as specified by ‘number of iterations for continuum point rejection’. For the initial fit and each subsequent iteration the terminal screen displays information about the fit. The right hand two columns give the number of channels included in the fit on each iteration, and the number rejected as being too far from the fitted curve. This number usually stabilises after a few iterations, i.e. the process converges to a limit and further iterations result in no further change. The first 8 columns of the display are the rms residuals of the (remaining) data points relative to polynomial fits of degrees 0 to 7. This full table is produced no matter what degree of polynomial fit you have selected. It enables you to see how low the degree of the polynomial can be without obtaining a noticeably worse fit to the data.
The continuum fit is drawn on the soft plot; if you have elected to write the continuum to an output Figaro file, the same data is written into the appropriate wavelength range of that file. It will overwrite any previous fit in this same ‘abline’ run which included the same wavelength range; i.e. all continuum fits in the one run go into the same output continuum file.
The next step is to use the cursor to specify the wavelength (channel) limits of the line. How this is done depends on the response you gave to the Limit parameter. The usual case is ‘limit=yes’: type ‘fit’ and centre the cursor on the last channel you want included, on each side of the line. Do the lower wavelength edge first.
The limits of the line are then drawn on the soft plot as vertical lines; for clarity they are drawn at the outer edges of the last channels included.
If the net area of the selected line is above the continuum rather than below it, the program reports that this is an emission line.
The program now has both the continuum and line specified, so it calculates and displays the quantities describing the line (for further details see Publ. Astron. Soc. Pacific 98, 1220 (1986)):
where F(k) is the data value in channel k and A(k) is the value of the fitted continuum there. The sum is over the specified channel range. A/ch indicates the dispersion in Angstroms per channel.
Define x(l) and x(h) such that
(For x(m), x(l), x(h) linear interpolation is used to find the fractional part of the ‘channel’ satisfying the desired equation.)
The width parameter is
For a line of Gaussian profile this gives exactly the FWHM, since x(l) and x(h) are at +/-1 standard deviation. For any reasonable centrally concentrated profile it gives a result very close to the FWHM. This definition is an attempt to obtain a parameter that is useful for a variety of line profiles. Note that the width parameter is considerably influenced by the cursor placement of the line boundaries, if these are at places where F(k) is not close to A(k). This is expected, and means that the width of the line as designated is partially under manual control. The line width parameter as defined above is less affected by low level wings far from the line centre than is the calculation of the width from the variance of the data about the mean, i.e. it refers more to the width of the core of the line.
The result is a dimensionless parameter (in percent), normalised by the line width. As an extreme case, a triangular line with a vertical low wavelength side has asymmetry of +8.1; if the high wavelength side is vertical, asymmetry -8.1. Higher values may result from deep narrow lines with prominent (noise) wings on one side. In general one expects the asymmetry parameter to be strongly affected by noise, and for many lines of low to moderate signal to noise it will be meaningless. Even the width parameter is rather noisy and should be interpreted carefully.
At this stage there is a short delay while the hard-copy plot is written to the designated file (see above for information about how to plot this after the run). The hard-copy gives a plot of the same wavelength range as selected and plotted on the soft plot; it shows the input data spectrum, the continuum fit, and two vertical lines bracketing the channel range selected for the line itself. A printed description on the same page as the plot lists all four parameters from analysis of the line (wavelength, equivalent width etc.) and also some of the details regarding the origin of the data and the type of continuum fit.
The program then loops back and asks for another centre wavelength for analysis. Give ‘quit’ to terminate if there are no other lines to analyse.
This section is more or less the documentation by Jeremy Walsh as included as on-line document in Figaro 3.0.
‘gauss’ is a Figaro program for interactive fitting of Gaussians to an emission or absorption line spectrum. The program is suitable for handling both low resolution data, where line fluxes of many lines are required, or high resolution data where a multiple Gaussian fit to a single line profile is obtained. The line-free continuum is fitted by a polynomial of desired order with iterative rejection of points greater than a given distance from the fit. The interactive line fit is optimised by minimising the sum of squared residuals between the observed and fitted profiles. The results of the fit are written to a data file. The fitting profile can be saved as a spectrum.
‘gauss’ requires an input spectrum in which the data points are equally spaced in wavelength. If you supply a spectrum with unequally spaced points it will abort with an error message. A spectrum with unequally spaced points can be converted to one with equally spaced points using Figaro command ‘scrunch’ (see Section 4.8.5).
Before running ‘gauss’ make sure that a soft and a hard plot device have been allocated using the Figaro commands ‘soft’ and ‘hard’. These are required for graphics (although use of the ‘hard’ device is not mandatory). The program will fail if ‘soft’ has not been specified.
The program will prompt for the following:
The spectrum is then plotted on the soft device in the lower box. If errors are available, then these are displayed in the upper box. The upper box is used for residuals on line and continuum fit. The scalar for the residuals indicates the ratio between the range of spectrum data to the range of residuals data, i.e. smaller errors and residuals have a larger scalar on the residuals box.
The continuum fitting menu next appears. The line or lines are demarcated from the continuum using the cursor, the parameters of the polynomial fitting set, and the fit performed. When satisfactory the fit and residuals are plotted ready for Gaussian fitting.
The options are as follows:
The Gaussian fitting menu now appears. The line or lines are simulated by one or more Gaussians.
The options are as follows:
For example: ‘P100’ shifts a Gaussian right by the profile extent; ‘H-20’ decreases the Gaussian height by 20% of the peak height of the original line; W10 widens the line by 10% of the total profile extent. At each change the old fit is erased and the new fit is redrawn.
S stops the alteration and re-plots the whole display, including residuals, on the latest fit. The rms (or mean fractional deviation in terms of the error bars) on the fit is written to the terminal.
The form of weighting for the residuals are prompted. Three options are available:
By employing weights more emphasis on points with either higher value or lower errors can help to constrain the fit. Method c is analogous to minimising chi-squared. Success cannot be guaranteed with the minimisation routine. In particular if constraints or chaining are employed then the optimisation can be more troublesome.
Usually a small change to the interactive fit will enable a minimum to be found if a failure occurs. However don’t expect miracles. The algorithm only finds a local minimum. If you don’t believe the fit, try another interactive fit and optimise again. The general rule is to try for a fit by a minimum number of Gaussians and to check that the deviations on the line fit are greater than or equal to the error bars, or the deviations on the continuum fit if errors not available.
If the optimisation is successful or the interactive fit is adopted, then position of peak, peak height, sigma and flux in each fitting Gaussian as well as the equivalent width on the whole profile and the simple (i.e. unweighted) rms deviation on the fit (mean deviation in terms of the error bars if errors are available) are written to the terminal and recorded on the data file (if applicable). On the ‘soft’ display, the fitting Gaussians, their sum and the residuals on the fit are displayed.
a) Low spectral resolution spectra (e.g. line strengths for single emission lines)
b) High spectral resolution spectra (e.g. radial velocity structure of a single profile)
The equivalent width tabulated is determined from the observed line profile and the continuum fit; it does not depend on any Gaussians fitted to the lines. The mean fractional deviation is defined as:
The section on échelle reduction is more or less the documentation “UCL Echelle Spectrograph Data Reduction Software” by William Lupton, dated 19 September 1988, as included as on-line document in Figaro 3.0. The appendices are not included here.
This note consists of a simple description of the current state of the échelle reduction software in Figaro. Wherever possible the échelle reduction package kindly provided by Jim McCarthy of Caltech has been used, but in several cases it has been necessary to modify the programs or else to use different approaches.
The software is now nearly in a state where it, together with standard Figaro applications, contains all the facilities necessary for the reduction of échelle data. Further development has been performed by the échelle data reduction programmer at UCL who produced a general purpose échelle data reduction package Echomop.
The diagrams illustrate the recommended processing steps starting from raw CCD or IPCS data. They are accompany a worked example.
Note that no special facilities are provided for bias subtraction, background estimation, sky subtraction or flat-fielding. It is assumed that you will use standard Figaro facilities for these operations.
Most of the rest of this note consists of an annotated example of how to proceed from a raw échelle image to a wavelength calibrated single merged spectrum.
Assume that you have just moved the spectrograph to a new configuration and that you have taken images of a continuum source such as a bright, smooth spectrum star (for tracking the échelle orders), of the Th-Ar lamp (for wavelength calibration) and of your object. In many cases the object spectra are sufficiently well exposed and free of discontinuities to allow direct order tracking. In such cases there will be no need for a separate continuum spectrum.
These files are in your current directory and are called:
In what follows it will be useful to refer to the diagrams in the overview section.
The first thing to do is to ensure that your data is in the conventional orientation with wavelength increasing from left to right and from bottom to top of the image and thus with order number decreasing from bottom to top of the image.
IPCS data will already be in this orientation but unfortunately CCD data must be rotated and flipped. This is quite time-consuming but is at present necessary for all CCD images. The following commands achieve this for the continuum:
with analogous commands for the arc and the object. The raw files can now be deleted if desired.
In what follows, assume that the IPCS raw data files ‘contraw.sdf’ have been renamed to ‘cont.sdf’ etc.
At present locating the orders is rather an interactive process (there is a program ‘echfind’ that does it automatically but it does not yet work very well in all cases and its use is not recommended).
First decide whether you are going to use the continuum source or the object to locate the orders. It is a good idea first to do a ‘ystract’ / ‘splot’ of the data to get a feel for the width, intensity and profiles of the orders. Sensible commands to use for the GEC chip and the continuum are:
Now display the image using ‘image’:
The next stage is to use ‘icur’ to define a point somewhere near the peak of each order that you wish to track and extract. You can track orders that are only partially on the image if you wish to but this is not recommended, since it could well affect the wavelength calibration, especially if only a small part of the free spectral range is being covered. It is quite important to choose points close to the peak intensity.
Now run ‘sdist’ to track the orders and fit polynomials to them:
The two non-obvious parameters are ‘trace’ and ‘width’. Specify ‘Gaussian’ for ‘trace’ if the profiles across the orders are roughly Gaussian and are not cut off by the dekker. Normally specify ‘COG’ otherwise but if there is a noticeable gradient along the profile you can try ‘Edge’ (they are identical in that both locate the rising and falling edges of the orders, but ‘COG’ estimates the centre by calculating the centre of gravity and ‘Edge’ estimates it simply by taking the mean of the edge positions). For ‘width’, specify an estimate of the FWHM for ‘Gaussian’ and an estimate of half the order width for ‘COG’ and ‘Edge’. If anything, underestimate it for ‘Gaussian’ and overestimate it for ‘COG’ and ‘Edge’.
Beware that sky data can confuse ‘sdist’ because it gives rise to profiles that don’t fit any of the trace modes. If this appears to be a problem, use ‘clip’ (which sets all data values below a given low value to that low value and sets all values above a given high value to that high value) to get rid of the sky data values that are causing the problem.
There is another program that may be useful here if you have moved to a new object that is not quite in the same place on the slit. ‘offdist’ operates on an ‘sdist.dat’ file and adjusts the constant terms so as to shift the tracked orders up or down by a specified amount.
You now have a choice. For quick-look extraction, you can create a mask image whose data values indicate which order (if any) each pixel belongs to. This mask is created by ‘echmask’ and applied by ‘maskext’. ‘echmask’ allows separate extraction of object and sky but requires the number of rows of object and sky data to be independent of order number. This is not acceptable when (as is usually the case) it is important to minimise the sky noise and to maximise the signal.
For final data reduction or where the use of a mask is not acceptable the orders can be straightened using ‘cdist’ (UCLES is extremely stable and preliminary results indicate that many images can be co-added prior to application of ‘cdist’, so the cost in processing time should be acceptable). Having straightened the orders, ‘echselect’ can be run to identify, for each order, the rows to be used for the object and those to be used for sky.
The result in both cases is a ‘collapsed échellogram’ where X is wavelength and Y is order number. Your object image will of course give rise to both an object and a sky échellogram.
As explained above, this method will not maximise signal to noise ratio, but it will do a reasonably good job, especially where the orders in question are not too far from being horizontally aligned on the detector.
‘sdist’ outputs an ‘sdist.dat’ file that contains details of the orders that it has tracked. This file is read by ‘echmask’, which produces a mask image that can be used for fast extraction of orders directly from (in the case of the CCD ‘irot90-ed’ and ‘irevy-d’) raw images. ‘echmask’ can cope with the case where the star / sky periscope is fitted and also allows you to specify the position of object and sky data relative to the centre of the order as determined by the tracking algorithm, but these details will be ignored here.
The normal straightforward behaviour is achieved by specifying ‘periscope’ false and giving zeroes for all widths and offsets. This causes a width derived from the ‘sdist’ fit to be used; if the fit was Gaussian the derived width is just twice the estimate that you gave to ‘sdist’ and if the fit was Edges it is the actual calculated width (the same width is used for each order and the third largest width of all the order widths is used so as to exclude atypical values). If you know better, you can specify your own value for ‘objwidth’—overestimate rather than underestimate so as to prevent noticeable jumps in the extracted data due to the slope and curvature of the orders.
The other thing that ‘echmask’ needs to know is the order numbers corresponding to the orders that it has tracked. The value that you give for ‘mstart’ is the order number corresponding to the first point that you selected with ‘icur’. There should always be an order near to the image centre. Don’t worry if you get it wrong—you can always adjust the order number by using ‘icadd’ to add or subtract the error from the mask structure. You have to add or subtract ten counts to adjust by one order, e.g. if the mask contains a data value of 420, this refers to order 42 and adding ten to it to make it 430 causes it to refer to order 43.
This is a typical run of ‘echmask’.
When running with the periscope, each order is split up into two parts, each of which looks rather like an order in its own right. If you are unsure how they are grouped, display an arc and all will be revealed. When using the periscope it is your responsibility to ensure that the first and second points selected with ‘icur’ correspond to the two parts of the same order and similarly for the third and fourth points etc.
The resulting mask image can be used for fast extraction of orders from (in the case of the CCD ‘irot90-ed’ and ‘irevy-d’) raw images taken at the same spectrograph configuration as it. This is done using the ‘maskext’ program. ‘maskext’ needs to be told the range of order numbers that you want to extract and this determines the Y size of the extracted file (referred to as a collapsed échellogram) and its Y units.
The ‘sub-order’ controls which bits of the order are extracted into the output image. A sub-order of 0 always extracts object and sky and sub-orders of 1 and 2 can be used to extract object and sky separately. If the periscope is fitted, sub-order 2 corresponds to the first encountered part of the order and sub-order 1 corresponds to the second encountered part of the order (so if the tracked orders went from the bottom upwards the data values in the mask are monotonically decreasing as you go from bottom to top). If the periscope is not fitted, sub-order 1 always refers to object and sub-order 2 always refers to sky. This is confusing and will probably change!
The trouble with curved orders is that when the object projects to only one or two pixels on the detector, the bulk of the signal will sometimes fall into one pixel and sometimes it will be split between several pixels. This means that there is no single correct number of rows of data to extract, forcing the extraction of unwanted sky as well as wanted signal.
Obviously it would be possible to provide a program that would make a sensible decision about how many rows of data to extract at each point along the order, and we intend investigating the use of some optimal weighted extraction scheme. However, any such scheme needs accurate knowledge of the noise characteristics of the data and it would take considerable effort to implement a reliable automatic extraction algorithm.
Accordingly, the recommended approach for accurate extraction is first to use the ‘cdist’ program to re-sample in the Y direction so as to straighten the orders. Experience shows that this program does an excellent job, with little discernible loss of resolution or variation of profile along the order. Once the orders are straightened, the number of rows to extract for object and for sky is merely a function of order number and not of wavelength.
The ‘cdist’ program uses the ‘sdist.dat’ file that was written by ‘sdist’. It uses the polynomial coefficients to re-sample in Y so as to straighten the corresponding orders.
Having got an image with straight orders, conceptually one wants to take a ‘ystract’ through somewhere near the centre of the orders, display with ‘splot’ and then, for each order, somehow identify which rows are to be used for object and for sky. Having done this, the relevant rows can be extracted into a collapsed échellogram of the same format as that produced by ‘maskext’.
‘echselect’ allows the user to indicate interactively the cross-sections of a corrected échellogram to be used as object and sky for the various orders. It then creates a collapsed échellogram for the object orders, and—optionally—one for the sky orders.
We are now at the stage where arc and object have both had their orders extracted into collapsed échellograms. These are two-dimensional images with X still being pixel number and Y being order number. The ‘echarc’ program is used for identifying lines in the arc. This writes wavelength calibration data as a two-dimensional arce.AXIS(1).MORE.FIGARO.DATA_ARRAY.DATA array and the arce.AXIS(1) structure is then copied to the obje.AXIS(1) structure using the ‘xcopy’ program.
‘echarc’ works by first of all doing the equivalent of the one-dimensional ‘arc’ program on a set of three or more orders that you nominate to be fitted interactively. Then it enters an automatic mode where lines are identified in all the other orders. It estimates wavelength in the other orders by fitting lines of constant (order number)*(wavelength) between the interactively fitted orders.
Experience with ‘echarc’ has shown that it is vital to include the extreme orders among those that are interactively fitted, that it may well be worth using four rather than three interactively fitted orders, that a good line list is absolutely vital and that if things start going wrong then they will probably stay wrong (use of Ctrl-C is recommended in this case). Care with the interactively fitted orders is usually rewarded.
The file ‘arlines.ech’ always contains details of identified lines and if a fit fails a good policy is to restart using the results from the previous fit and perhaps selecting slightly different interactively fitted orders, e.g. add the orders that were worst in the automatic phase from the previous run. All the lines previously identified in the orders that are to be interactively fitted are still available so a re-run is not too time-consuming.
Also note that it is necessary to rename any old files of name ‘arlines.ech’ before running ‘echarc’.
The ‘dowaves’ keyword should always be false (if it is true then the wavelength information is written to a separate file and we want it in the input file for the purpose of scrunching). Also note that if the ‘monitor’ keyword is true then a graphical record of how the automatic mode is proceeding will be output to the soft device.
Here is a typical example of ‘echarc’.
After a successful arc fit, use ‘xcopy’ to copy the ‘arce.AXIS(1)’ structure to the ‘obje.AXIS(1)’ structure.
The next step is to scrunch the wavelength calibrated ‘obje’ and ‘arce’ files. Contrary to possible expectation the program used is ‘scrunch’ rather than ‘iscrunch’. Rather than using an ‘arlines.iar’ file, ‘scrunch’ uses the two-dimensional .AXIS(1).MORE.FIGARO.DATA_ARRAY.DATA array produced by ‘echarc’ as its source of wavelength information. You tell it the wavelength range and the number of bins that you want to scrunch into and whether to use a linear or logarithmic wavelength scale and it does the rest.
Since the échelle spectra are such that a pixel corresponds to a fixed velocity interval irrespective of wavelength (e.g., 22 microns = 2.73 km/s), a logarithmic wavelength scale will give bins that each correspond to the same number of pixels, whereas a linear wavelength scale will give bins at the red end that significantly oversample the resolution element relative to bins at the blue end. However, a linear scrunch will still probably be the favoured option in many cases.
It’s probably worth using ‘hdstrace’ to look at the file to find out the minimum and maximum wavelengths before running ‘scrunch’.
Every order will be scrunched into the full wavelength range (4550 to 5000 Angstroms in the above example) and consequently most pixels will be zero.
The final stage is to merge the scrunched orders from possibly several images into a single long spectrum. This is done using the ‘echmerge’ program. If the data is from a region of the field where the full free spectral range is obtained then adjacent orders will overlap and the contributions from the overlapping orders are weighted with estimates of their inverse variances on the assumption of purely Poisson noise. Where one contribution is below a given fraction of the other then it is ignored completely. The variance estimates are based on median filtered versions of the input orders and the median filter box size is controlled by the ‘box’ parameter.
The input images must all have the same X size, units and values and can be one-dimensional or two-dimensional. The output image is always one-dimensional and can be the same as either of the input files, in which case no new image is created. The second input image can be specified as blank in which case it is not required or used.
Here there is just a reference list of the programs that were run above.
3J.B. Oke and J.E. Gunn, 1983, Astophys. J, 266, pp713-717.
4J.B. Oke and J.E. Gunn, 1983, Astophys. J, 266, pp713-717.
5A.V. Filippenko and J.L. Greenstein, 1984, Publ. Astron. Soc. Pacific, 96, pp530-536.
