This section contains articles written by Rachael Padman for the JCMT Newsletter and are included here since they may be difficult to obtain from other sources.
Traditionally, radio spectral-line analysis programs have displayed velocities with respect to the Local Standard of Rest (LSR), using the radio definition of velocity. It is now possible to observe at JCMT using Heliocentric, Barycentric and Telluric velocity frames, and optical and relativistic velocity definitions, and SPECX V6.2 (now released to Starlink) contains the appropriate facilities for dealing with such data.
Before describing what SPECX does, it is necessary to review, briefly, what is required (see also Gordon 1976). We all know that any relative velocity between source and observer gives rise to a corresponding Doppler shift in the received frequency. For small velocities, then, radioastronomers (who think in frequencies) write
| (1) |
where is the velocity away from the observer. Likewise, optical astronomers write a similar formula in wavelengths:
| (2) |
Finally, it seems that extragalactic observers believe in special relativity, which of course tells us that:
| (3) |
(ignoring the smaller doppler shift due to any transverse component of velocity). By a simple application of the binomial theorem, we can see that these are all equivalent for . A problem arises because these equations are used in reverse to define a velocity, and this definition is often used where the approximations in (1) or (2) are inadmissible. Even the (special) relativistic formula applies only in local frames, which are not the same thing at all as distant cosmological frames. So none of these formula are “right”.
Normally we reduce our velocities to one of a few commonly agreed standards: Radio astronomers prefer the Local Standard of Rest (LSR); extragalactic astronomers apparently prefer a Heliocentric system, whilst for some purposes, such as instrument or atmospheric diagnostics it may be the Telluric frame that is of most interest. Or indeed we may have some favourite other frame which is defined with respect to one of these "standard" frames. One that is used quite often is the source frame – for example, if we want the frequency scale to be that appropriate to a non-moving source we have to transform all telluric (observed) frequencies to a frame comoving with the source. In the case of Orion for example this is 9 km/s wrt the LSR, and few of us know the heliocentric velocity.
Thus our velocity scale is specified by four parameters:
Here I use the general term "frame" to imply this full specification.
There are two separate issues. First, we need to observe such that the line of interest actually occurs within the passband of the instrument (normally, but not always, at the centre of the passband). For a heterodyne receiver we do this by adjusting the local oscillator frequency by an amount equal to the total doppler frequency shift expected for the line. Second, we then need to display the spectrum. Unfortunately, the spectrometer is back on Earth, at the telescope, so we need to map the I.F. frequencies of individual spectrometer channels onto velocity (or frequency) space in some other chosen frame. This ‘display’ frame need not necessarily be the same as the ‘observation’ frame. For example, whatever the observation frame, we may wish to display the output in one of several standard frames:
The solution adopted in SPECX is to view all calculations of velocity and/or frequency as a two stage process. In the first, the spectral header information is used to calculate the telluric centre frequency of the observation. That is, we deduce the true frequency as measured at the telescope of a signal appearing in the centre channel of the spectrometer. SPECX then produces an ‘X-array’ which contains, for each spectrometer channel, the telluric frequency corresponding to that channel. Finally, this array is transformed back to the display frame. By default the display frame is the same as the observation frame, which is encoded in the SPECX and GSD scan headers, but it may be changed to any other frame if you like.
To summarize: SPECX will normally display the data using the velocity frame in which it was observed; however you can use SET-VELOCITY-FRAME to select another frame, and optionally you can use SET-LINE-REST-FREQ to choose another reference frequency for the velocity transformation. Current options for frames include LSR, Geocentric, Heliocentric and Telluric, and you have a choice of Radio, Optical and Relativistic velocity laws. Because of the bewildering number of combinations of these variables, the header on the X-axis of the plot has been modified to give a full specification.
The command CHANGE-SIDEBAND does the following: First, it calculates the telluric centre
frequency. Then it adds or subtracts twice the I.F. (depending on whether the observations are in LSB
or USB respectively). Finally, it transforms this telluric frequency to the nominal frame of the
observation, files it as F_CEN in the scan header, and changes the sense of the frequency
step in the spectrometer (F_INC). GSD spectra stored with storage task V6 or earlier do
not have the local oscillator frequency stored, so SPECX instead asks you for the current
sideband and I.F.; for V7 data or later SPECX can deduce these from the data in the GSD scan
header.
Let us take a specific example. We have a RxB3i spectrum which is centred on CS 7-6 (rest freq. 342.883 GHz) in the lower sideband. The nominal I.F. is 1.5 GHz. In order to look at absolute "rest frame" frequencies, we set the display frame to be that of the source itself:
(see Fig 1). (There is a predefined variable "vlsr" equated to the actual value for the current spectrum, so we can quote this in response to the last query rather than looking up the actual value.)
For an AOSC spectrum, we would now regrid to a linear scale (see next section for more details on this). Then we go ahead and change the sideband:
As suggested, change the reference frequency to CO 3-2:
If we have done all this correctly we now get the absolute frequency plot shown below
(Fig 2). A macro has been provided to do all this – just type IMAGE-FREQ (this calls
specx_command:image.spx).
It is well known that the frequency scales of Acousto-Optical spectrometers tend to be non-linear. To
this end, the SPECX frequency calculation stuff contains a facility for correcting the frequency using a
polynomial fit to the frequency error. For AOSC, the cubic term dominates, and the maximum error is
about 0.6 MHz, but for somewhat complicated reasons, there is an additional zero-order
(d.c.) term of some 11 or 12 channels, or about 3 MHz. The implication of this is that a line
you might expect to come out in the centre channel (1024.5) of the spectrum will actually
emerge some 3 MHz away, and to compensate for this the GSD header files correct the
reference-frame observed frequency, F_CEN, appropriately. This makes the display work fine, but
has implications for when we want to look at image frequencies (as will be shown in a
moment).
To correct the frequency scale for AOSC data in particular, I have written a small macro, FRQFIX.SPX, which is kept in the "standard" command-files directory. It is invoked by the symbol LINEARIZE-AOSC-FREQ:
This macro in fact also removes the 3-MHz offset from F_CEN, instead adding it into the zero-order
term of the correction polynomial, which simplifies calculation of the I.F. LINEARIZE also does
various checks and tidies up the program flags. If you choose to REGRID the data at this point you
transform the data to a truly linear scale (in current units). Otherwise the data will still
display correctly, as long as you have the linearization turned on (set fcorrect=true or use
SET-X).
If you now want to look at the other sideband of AOSC data, remember that the effective I.F. at the
centre of the passband is not what you might think it is, but is actually 1.503 GHz (or 3.943
GHz for RxA1 or RxB2). If on the other hand you apply the LINEARIZE-AOSC macro
the data will be regridded onto a linear scale and shifted to the nominal channel. In this case
the I.F. is the standard 1.5 or 3.94 GHz. There is one Awful Warning: you cannot apply
linearization to data after you have done a CHANGE-SIDE. That flips the spectrometer frequency
scale, and so will apply the correction ’the wrong way round’. So after you have done a
CHANGE-SIDE, make sure that linearization is turned off, either through SET-X or by typing
"fcorrect=false".
I thank Per Friberg, Goeran Sandell and Chris Mayer for their help in untangling this mess, Louis Noreau for pointing out the need for other velocity frames and scaling laws, and Paul Feldman for bringing other infelicities in the velocity and frequency scaling to my attention. The system in SPECX V6.2 has been written from scratch, and so does not have 10 years of testing behind it; please bring any demonstrable faults to my attention. (I do not have Donald Knuth’s confidence, and do not offer a steadily increasing reward for each bug found in my code.)
Finally, a brief note about map-making. People are making some really quite big maps. SPECX tries to hold the entire cube in virtual memory at one time, and in fact with various permutations of INTERPOLATE and ROTATE, may want to have up to 3 copies of the cube resident. If excessive paging is to be avoided, then it is desirable that all these fit into the available physical memory of the machine (the working set). The usual consequence if they do not is that the whole machine grinds to a halt, without any apparent error...
You can minimize the amount of virtual memory required by only mapping those spectral channels of interest. Use TRUNCATE (or DROP-CHANNELS) to dispense with spectral channels lying far away from the line; use BIN-SPECTRUM where you have higher resolution than you need. Even if you have plenty of memory to spare, you can greatly speed up the map-making process by using as few channels as possible. To reinforce this point, when you do an OPEN-MAP you are now asked explicitly for the number of spectral channels in the cube. You can also eliminate one cube from memory altogether by using ‘interpolation-on-demand’ (an option in INTERPOLATE-MAP), albeit at the cost of slightly increased time to make any particular map.
As noted in the manual, SPECX maps are really channel maps, and no information is stored in the header about individual spectra. Thus if you want to map data taken in the other sideband from that of the map header, or taken with a different spectrometer etc, you first need to INVERT, SHIFT, TRUNCATE etc to get your spectrum into the right form. There is a new macro that does this automatically: it can be invoked by the command symbol CONVERT-TO-MAP-FORMAT this macro is stored along with all the other standard macros in the directory with logical name "specx_command".) Just use this command before you ADD-TO-MAP to convert your current spectrum to the same frequency scaling etc as the map header.
M.A. Gordon. Chap 6.1 in "Methods of Experimental Physics", vol 12-C, (Astrophysics – Radio Observations). Academic Press, NY, 1976
Following publication of the last issue of the Newsletter, Pat Wallace (author of the SLALIB package) wrote to point out that my use of the phrase “the LSR” may have been confusing. The “local standard of rest” is strictly a defined frame, but is meant approximately to represent the basic solar motion with respect to the neighbouring stars. This motion is established by measuring the radial velocities of the stars in the solar neighbourhood. The number you get depends on the depth of the sample (how far out you go), and the spectral classes of the stars you use. Thus there are several determinations of the kinematic solar motion, which differ by a few kilometres per second, and by a few degrees in direction. Spectral line radio astronomers however traditionally use the defined “LSR”, which has a velocity of 20km/s in the direction of 18h,+30d(B1900), and differs slightly from the most modern determinations of the solar motion.
Pat has recently altered the SLALIB routine SLA_RVLSR to provide a better estimate of the
dynamical solar motion — i.e., the motion of the Sun with respect to the appropriate
circular orbit around the galactic centre. This also makes SLA_RVLSR consistent with
the routine SLA_RVGALC, but unfortunately there is now a velocity difference of up to
km/s
between the velocities calculated by SLA_RVLSR and those used by radio astronomers (and SPECX). I
propose not to implement this new definition in SPECX, as I suspect it will only cause confusion if
Orion starts to come out with "an lsr velocity" of 6km/s. Pat Wallace’s program RV, which is
included in the JCMT utilities, is based on the SLALIB routines, but the older version is
being retained to prevent undue confusion arising as a result of the change in the SLALIB
routines.
One further note: SPECX and the JCMT control system have both in fact been using a velocity of 20km/s towards an apex of 18h,+30d(B1950), as used in the original Bonn software. For consistency I will change the epoch used in SPECX to B1900 in the next release, but the difference in velocities will be very small («0.1km/s), so shouldn’t be observable for AOSC data. The apex used in the telescope control software will also be changed to B1900 as soon as practicable.
My thanks to Pat Wallace, Chris Mayer and Per Friberg for helping me to sort out what was actually going on here.
One serious bug has surfaced in SPECX6.2, which means that some maps are “unopenable”. The program starts to open them, and then falls over. A kluge which fixes this in some cases is to open another map first, then open the offending map afterwards. However the only real fix is to replace the incorrect routine. Any site which has problems can obtain a copy of the routine direct from me, or from John Lightfoot at ROE (REVAD::JFL).
A slightly revised version of SPECX (V6.2-A) is running at JCMT. This fixes the bug described above, and one or two other very minor problems. Changes include now being able to open up to 8 SPECX data files (was 3), and being able to change the parameters controlling the “MRAO colour spiral” (colour 5 for colour-greyscale mapping). There is also a tentative implementation of logarithmic greyscales – just hit the ‘0’ key to toggle from log to linear and back.
I have added a command definition UTILS, which prints a file that describes the “system
macros” – .spx files written to accomplish certain oft-needed functions, including those
such as ‘map-average’ and ‘fetch’ that have been described in previous versions of the
newsletter.
The good news (for those for whom that sort of news is good) is that, thanks largely to the efforts of Horst Meyerdierks at The University of Edinburgh, and following an intensive burst of work over Christmas on my part, SPECX should shortly be available on Starlink platforms. Large parts of the port are complete but it is not quite clear when there will be enough of a system to be worth releasing it to Starlink. At the date of writing (4th January) neither GSD nor FITS i/o has been implemented, so there is still no way to import data other than by direct conversion of VAX files in SPECX internal format. Horst has written a standard SPECX command to do this - for the present anyway it produces .SDF files, which use Starlink NDF and HDS libraries. Remo Tilanus has been working on a GSD reader, and converting the FITS stuff should not be too difficult (although tape i/o will probably not be supported under ), so with any luck it will not be too difficult to finish off V6.4.
Unfortunately a number of the nicest features of the VAX version - such as CTRL C handling and proper support for dual-screen alpha graphics terminals - are unlikely to be available in the first version (6.4), and users may notice some "debug" type statements which are unfortunately necessary to circumvent bugs in the f77 compiler, but for most purposes SPECX V6.4 should satisfy the need. With any luck it will be released to Starlink during February, and if not, then soon after. Whether the documentation will be available by then is another matter.
Horst is aware of the need to be able to export SPECX to non-Starlink sites without them needing to take on the whole of the Starlink environment, and I am assured that this will be possible. However in the first instance it is probable that software will be available only through ROE, or later through the Starlink Librarian — in any case, probably not from me directly.
Meanwhile, for those who do have access to a VAX still, and while you are at the telescope, the improvements in V6.3 should make life easier. In fact the release of this version has been prompted mainly by the successful commissioning of the DAS, which in turn has required major changes in the internals of SPECX to accommodate extra header information required to define the frequency scales accurately. Several new commands have been added for dealing with DAS data, while others have been modified slightly.
Most importantly, the SPECX internal datafile and map format has changed. As always, SPECX 6.3 is able to read any existing files and maps you may have, but it is no longer able to write to them. Map files (with the .MAP extension) are converted to V3 the first time they are opened by the new program — if you intend also to use earlier versions of SPECX then make a copy of the map before you use V6.3 for the first time. This is not true of ordinary datafiles; you just can’t write spectra to them any more, but it is a simple matter to open a new file and then copy the spectra across one by one or in a DO loop.
The reason these file formats have had to change is to accommodate the LO and IF frequencies in the header, along with the velocity of the source in various frames. These were not previously available in the GSD file, but were derived within SPECX from the specified source-frame centre frequency of each sub-band. This is not however a natural way to do things, and proved prone to error with increasingly complex data such as those produced by the DAS. Hence the change. Since things were changing I have taken the opportunity to tidy up some other relics of bygone computer systems — RA and DEC are now stored internally as R*8 degrees, and the map offsets are stored as reals. This has meant changes to some of the macros supplied with SPECX, and if you use these variables in your own macros then you will need to modify them accordingly.
Bad-channel handling has been available for some time in the various mapping functions, although it tends to go unseen (except when you forget to do an INTERPOLATE for sparsely sampled data). A version of this has now been introduced for single-spectrum data. Channels in spectra set to this value are not displayed on plots, or used in computing quantities such as maximum intensity or line width. Channels which end up being unspecified, as the result of a SHIFT operation say, are now set to the bad channel value, rather than zero as before. Bad channels can be interpolated over with smoothing commands such as HANN, SMOOTH, CONVOLVE, BIN etc.
You can determine the current setting of the bad channel value by doing:
and if you don’t like it, you can change it by doing:
for example. The one place where bad channels in spectra are not handled well is in spectra in a map
cube — fewest problems will result if you smooth over such bad channels before ADDing to the map,
or choose a value of badpix_value close to zero (I tend to use 1.e-5 by default). Although we did have
some teething problems with this feature, it seems to be working pretty well now — It is particularly
useful for dealing with DAS data, where the end channels of each subband tend to be set to
or
some such similar value, and totally destroy the auto-scaling of plots.
Note that if you change the value of badpix_value as shown above, then previously “hidden”
channels in your data will now be displayed, and/or used in reduction operations.
A few mostly cosmetic changes have been introduced here. Displays involving greyscales will mostly now have a scale-bar displayed wherever there seems to be most room. This is not foolproof however, and for some unfortunately sized maps the scale-bar may vanish off the edge of the plot area – just do SET-MAP-SIZE if this is the case. Note that the scale bar will not appear in the case of PLOT-LINE-PARAMETERS when more than one map is made — the various maps in this case have different scales, and it did not seem sensible to try to generate scale-bars for each independent panel.
Also, you can now display RA and Dec coordinates (only) in all forms of map display (GRID, CHANNEL-MAP, GREY, CONTOUR and PLOT-LINE-PAR) in standard sexagesimal notation. That is, absolute positions are displayed, in HMS and DMS respectively. You can however still display things in the old way if you want – the map centre RA and DEC are now indicated in the axis labels. Use SET-MAP-PARAMETERS to change the default axis scaling.
There have been ongoing problems with exchanging single spectra in FITS format between SPECX and CLASS. We (Remo Tilanus and I) believe that this has now been cured, following minor changes to both programs (with the cooperation of the CLASS authors at IRAM). Certainly SPECX now seems to produce valid CLASS spectra. I have also resurrected John Richer’s FITS reader, and with a great deal of hacking turned it into a standard command that reads disk-FITS spectra (SIMPLE=.TRUE. only) into the X position of the stack. This command is READ-FITS-SPECTRUM — open the FITS file with OPEN-FITS-FILE and close it afterwards with CLOSE-FITS-FILE as before (note that the names of these commands have been changed from OPEN-FITS-OUTPUT-FILE and CLOSE-FITS-OUTPUT-FILE respectively).
FITS maps are more of a problem. Claire Chandler has recently modified WRITE-FITS-MAP to make it compatible with AIPS (it always used to be, so I assume that AIPS has changed too...), but this is probably not compatible with CLASS. We may eventually need separate commands to write FITS maps for different targets.
It is now possible to produce plots which have the image sideband frequency along the top axis. To get this, all you have to do is to choose the absolute frequency options in SET-X-SCALE. Note however that it is unlikely to be correct unless the display frame is the one in which the source is at rest. Since you *may* have observed with (say) the LSR velocity not equal to that of the source itself, I can’t protect you from this. But in most cases it will be sufficient to do
That is, choose a velocity reference frame which is defined as having velocity offset with respect to the standard of rest the same as that of the source itself (modify as appropriate for heliocentric, geocentric or telluric frames, and for different velocities or velocity definitions).
HOWEVER... You must realize that most plotting packages, PGPLOT amongst them, accept REAL*4 numbers only. So if you display things in absolute frequency, it is not possible to get velocities accurate to better than a couple of MHz in three hundred odd GHz, no matter how accurately they are calculated internally in SPECX. For more accurate displays, choose a convenient reference frequency in the absolute frequency option in SET-X. For example, if you have a set of lines in the range 338 – 339 GHz, then choose a reference frequency of 338 GHz, and remember to add this back on to any frequency you determine with FIT-GAUSSIAN, or using the cursor.
My thanks to Remo Tilanus in particular for his help, provision of DAS data, and rapid turn around of queries while we were solving the problems associated with the DAS. My apologies to anyone who may have been inconvenienced by problems in the Beta-test version — I hope you appreciate that it was the only way to make DAS data available at all.