 News |
FAQ
|
| Search
|
Home
|
| Getting Started | Documentation | Glish | Learn More | Programming | Contact Us |
|
| Version 1.6 Build 363 |
|
Dana S. Balser, J. A. Braatz, & Joseph. P. McMullin
This chapter describes how to reduce and analyze GBT spectral line data using AIPS++. Although the examples are tailored for the GBT it should be relatively straight forward to modify them for any telescope. This chapter is expected to evolve with time as users begin making more observations with the GBT and provide feedback on the data reduction software. Comments are appreciated.
In §9.2 the data flow from the GBT to AIPS++ is described. In some circumstances the raw data require examination and flagging before calibration. Tools for data examination and inspection are discussed in §9.3. Currently GBT data are not calibrated automatically. This is expected to change soon. In §9.4 several user Glish functions are provided which are used to produce a crude calibration. The intention of this section is not to discuss the details of calibration but to provide the user with examples on how data can be processed in AIPS++ at the Glish level. Some basic data analysis examples are given in §9.5. Conclusions are in §9.6.
The data used in the examples below were taken from the GBT L-band (1.15-1.73 GHz) receiver when the telescope was fixed in (Az, El); the contractors were still working on the telescope. The spectral processor (SP) was used as the backend.
The GBT monitor and control system coordinates an observation and writes data from various devices into several FITS files. A filler combines the information in the FITS files into what is called a measurement set (MS). This operation can be performed off-line for the GBT project AST_SPEC_BASEL4 using the following shell command:
gbtmsfiller project=/home/gbtdata/AST_SPEC_BASEL4The resulting MS is called AST_SPEC_BASEL4_SP by default. The ``SP'' notes that the backend used was the spectral processor. The MS consist of a collection of AIPS++ tables and subtables.
To view the MS start an AIPS++ session. Here we use the stable version of AIPS++. The first line executes a script which sets various paths, environment variables, etc.. The second line starts a session of AIPS++. A ``Log Messages'' and a ``Tool Manager'' window should appear. In the window in which you started AIPS++ you should see the Glish prompt ``-''. For information about the Log Messages and Tool Manager windows see ``Getting Started in AIPS++''9.1
source /aips++/stable/aipsinit.csh
aips++
reading .glishrc_local from /home/aips++/stable
************************* Welcome to AIPS++ *****************************
Copyright (C) 1995-2000,2001 Associated Universities, Inc. Washington DC, USA
AIPS++ comes with ABSOLUTELY NO WARRANTY
Version 1.5 (build 284) on solaris_egcs
for more details, type about()
Loaded system packages: utility
Type help() for help
Glish version 2.7.
-
The MS can be viewed by using the table browser. First set an arbitrary variable (e.g., ``t'') equivalent to the table and then use the browse function. A table window should appear with columns and rows. In general there are additional subtables which can be viewed by selecting ``table summary'' under the ``view'' menu.
- t := table('AST_SPEC_BASEL4');
- t.browse();
The observations were made with the spectral processor (SP) configured in the 2 x 1024 mode. That is, there were two signals transmitted from the receiver to the backend with 1024 channels each. These signals have been mixed to intermediate frequencies called IFs. In this case the first IF corresponds to a right circularly polarized (RCP) signal and the second IF to a left circularly polarized (LCP) signal. The bandwidth for each IF was 40 MHz. The RCP signal was centered on 1410 MHz, while the LCP signal was centered on 1380 MHz. The SP integration time was set to 60 s with a duration (or scan time) of 10 min.
Currently the GBT filler does not calibrate the data. So the intensity scale is in counts which are proportional to the total power. A calibration noise diode with a measured intensity was fired ON and OFF with a cycle period of 1 s. Therefore there are two separate phases per IF. These phases have been averaged separately over each integration period. Therefore, each scan consists of 2 IF x 2 phases x 10 subscans = 40 spectra or records. The entire MS consist of scans from 181 to 240 for a total of 60 scans or 2400 records. To summarize:
Record Scan Subscan IF Phase 1 181 1 1 on 2 181 2 2 on 3 181 3 1 off 4 181 4 2 off 5 181 5 1 on 6 181 6 2 on 7 181 7 1 off 8 181 8 2 off . . . 33 181 33 1 on 34 181 34 2 on 35 181 35 1 off 36 181 36 2 off 37 181 37 1 on 38 181 38 2 on 39 181 39 1 off 40 181 40 2 off 41 182 1 1 on . . .
In the future a flag will exist in the filler that will allow the user to calibrate the data and produce antenna temperatures based on a calibration table. However, there may be some instances when processing the raw data is desirable. For example, it may be best to flag RFI before the calibration is applied.
The data can be examined in various forms and flagged using the AIPS++ tool msplot. Include the tool and then create a version of msplot using the MS as input.
- include 'msplot.g';
- mp := msplot('AST_SPEC_BASEL4');
A msplot window should appear. Select the appropriate data subset
using the selection criteria. Notice that line plots as well as
images can be viewed. Use the cursor to flag regions of an image,
if desired. The actual data are not removed but a flag table
which acts on the data is generated and can be applied to the data.
This section describes some simple Glish functions which were used to calibrate the data. The main purpose here is not to discuss how GBT data are calibrated but to illustrate how AIPS++ can be used to process a MS. This will be done primarily using the scan-based single-dish analysis tool (SSA). This tool is currently under development and will evolve with time.
The SSA is a tool which is a layer on the single-dish package called DISH (see NOTE 236 - DISH: User's Manual).9.2 DISH can be accessed either at the graphical user interface or at the command line interface. The SSA was developed to provide simple, compact functionality at the command line.
First this tool needs to be included in your AIPS++ session. Next create an SSA session. Note that the variable ``gbssa'' is arbitrary and just stands for Green Bank's SSA. The main DISH window should popup along with the DISH plotter.
- include 'ssa.g' - gbssa := ssa();
Currently there is no formal documentation for the SSA. Nevertheless, the functionality within SSA can be listed by using the Glish function field_names.
- field_names(gbssa) type gui average baseline clearrm dish filein fileout gauss getscan help history listscans listsummary open plotscan plotxy rmadd save scansub scandiv select scale smooth statistics writeA quick help for each function can be listed by using the help function:
- gbssa.help(F);
## ssa = scan-based single dish analysis
## This is a compacted CLI syntax for more facile scripting and
## command line operation. Motivated by GBT observers.
##
.
.
.
##write Description: Write X,Y axes to an ascii file on disk.
## Example: myssa.write('./testwrite');
## Returns: F
## Note: It comments:
## Spectrum written to file: ./testwrite
## Produces: A file on disk, in the current directory called
## 'testwrite'.
T
To list the code for each SSA function just execute the function
without any parameters or parentheses. For example, to list
the Glish code for the open function:
- gbssa.open
function (val file, val readonly = T, val new = F) {
{
{
{
{
(ok := sd.open(file, readonly, new = new));
(size := sd.rm().size());}
(myws := sd.rm().getnames(size));}
(private.wsname := myws);}
(ok := public.filein(private.wsname));}
if is_fail(ok) {{ return throw(FAIL: Could not open file) }}
else {{ return T }}}
-
In order to work on a MS we have to load it into DISH. For historical reasons DISH does not currently work with a measurement set (MS). Instead the MS is converted into a single-dish (SD) record or working set. The MS is a general data format for radio astronomy which includes aperture synthesis as well as single-dish data. The SD record is specifically designed for single-dish data.
The SSA open function will perform this task. Since many working sets may be loaded into DISH another SSA function filein is necessary to point to the appropriate data. Alternatively, the MS could have been opened by selecting ``open'' under the file menu in the DISH window and then using the cursor to select it.
- gbssa.open('AST_SPEC_BASEL4')
- gbssa.filein('AST_SPEC_BASEL41');
Next are some example user scripts enclosed in boxes which use the SSA and basic Glish functionality. These scripts are included just like any other Glish file.
The first function is called tsys which calculates the system temperature. The scan and subscan numbers are required, along with the IF channel, the total number of phases, and the temperature of the noise diode in Kelvin. First the noise tube state for ON and OFF are selected using the SSA getscan function; then the system temperature is calculated for each spectral channel; lastly the mean system temperature across the band is returned.
- tsys(181, 1) 15.4910584
The function tant is similar to tsys except the whole array is returned instead of the mean. Actually, the entire SD record is returned and contains header information as well as the data array. Notice that initially the variable tant is set equivalent to the on scan. This creates all of the SD record fields. Then only the data array is redefined.
- scan181_1 := tant(181, 1) - field_names(scan181_1) data header hist otherIf we want to access the source name, which resides in the header field,
- scan181_1.header.source_name testTo print the entire contents of the header:
- scan181_1.header [time=[type=epoch, refer=UTC, m0=[value=51935.8791, unit=d]], scan_number=181, source_name=test, direction=[type=direction, refer=J2000, m1=[value=1.57079633, unit=rad], m0=[unit=rad, value=0]], refdirection=[type=direction, refer=J2000, m1=[value=1.57079633, unit=rad], m0=[unit=rad, value=0]], veldef=RADIO, transition=, exposure=30.149786, duration=59.9992332, observer=, project=, resolution=1000000, bandwidth=40000000, tcal=0, trx=0, tsys=0, telescope=GBT, telescope_position=[type=position, refer=WGS84, m2=[value=807.435, unit=m], m1=[unit=rad, value=0.670784476], m0=[unit=rad, value=-1.39346797]], pressure=0, tambient=0, dewpoint=0, wind_dir=0, wind_speed=0, azel=[type=direction, refer=AZEL, m1=[value=0.667359901, unit=rad], m0=[unit=rad, value=-5.0113769e-05]]] -
The function calout calibrates the data and saves the results as a new data set. This function is essentially a big loop which calculates the antenna temperature using tant and then uses the the SSA function save to write data to a new working set (WS).
In general using big loops is not good programming practice in Glish since their performance is poor. Before calout can be used a new working set (let's call it spdata) must be opened. Notice that we have to get permission to write to the working set and specify that it is new. The default is readonly and old.
- gbssa.open('spdata', readonly=F, new=T);
- gbssa.fileout('spdata');
- calout(181, 240);
181 1 1
181 1 2
181 2 1
181 2 2
.
.
.
240 9 1
240 9 2
240 10 1
240 10 2
T
In this section several very basic data reduction procedures are demonstrated. First re-open the working set ``spdata'' with readonly protection, the default.
- gbssa.open('spdata');
- gbssa.filein('spdata')
A summary of the scans in the working set (WS) can be listed by using
the SSA function listsummary. Alternatively, the browse button on the
DISH window can be used to popup another window with a summary of the
records as well.
- gbssa.listsummary(); 181 test 26-Jan-2001/21:05:56.152 181 test 26-Jan-2001/21:05:56.152 . . . 240 test 27-Jan-2001/08:13:17.015 240 test 27-Jan-2001/08:13:17.015 T
Notice that now there are a total of 20 subscans for each scan number, for a total of 1200 records. A given record can be plotted by either clicking the specified record from the DISH browser or by using the SSA function plotscan.
- scan190_8 := gbssa.getscan(190, 8); - gbssa.plotscan(scan190_8); T
A common task in any data reduction package is to average spectra. This can be performed with the graphical user interface (GUI) in DISH by selecting the appropriate item under the operations menu. In this case both ``selection'' and ``averaging'' are required. Notice that as each operation is selected that the GUI expands to include the new functionality. First the appropriate spectra must be selected and then averaged using some criteria. For example, how should the spectra be aligned or weighted for the averaging process? To perform the average just click on apply.
Alternatively, a group of records can be averaged together by using the SSA function average. For example, let us average all of the subscans for scan 190, IF 1. Recall that the IF's alternate so the odd subscans correspond to IF 1 and the even subscans correspond to IF 2. The subscans can either be specified directly or by using the Glish function seq. After the average plot the results.
- scan190 := gbssa.average(190, [1, 3, 5, 7, 9, 11, 13, 15, 17, 19]); - scan190 := gbssa.average(190, seq(1, 20, 2)); - gbssa.plotscan(scan190);
In practice the user may want to inspect each spectrum before the averaging process. The user function avgscan allows the user to inspect each spectrum before it gets placed into the group to average. Here each IF is done separately.
The following demonstrates how to run the function avgscan to view and average only the first 4 subscans of scan 190. The ``> >'' prompt will wait for user input. Just hit return if the displayed spectrum is okay to average and ``n'' if not. The resulting averaged spectrum will be displayed.
- avgscan(190, 190, 4) Average IF 1 subscans for scan 190 scan 190 1 >> scan 190 3 >> Average IF 2 subscans for scan 190 scan 190 2 >> scan 190 4 >> T
A common data reduction procedure in spectroscopy is to model the spectral baseline. The spectral line might lie on a much larger continuum background which must be subtracted to measure the line intensity. In many cases the underlying continuum shape is modified by instrumental effects, such as standing waves caused by reflections from various parts of the telescope structure.
From the DISH GUI select baselines under the operations menu. The type of fit (polynomial or sinusoid), the order of the fit, and how to apply the fit (show or subtract) can all be selected. The spectral regions or ranges to model can be defined by typing them directly in the field box or by using the cursor.
The user function basefit will plot a specified scan, prompt the user for the line free regions to be fit by a polynomial, calculate and plot the polynomial model, and lastly subtract the model from the data and plot the residuals.
The following example will fit a 3rd order polynomial to scan 191, subscan 1. Notice that the baseline region ranges are defined as a series of Glish records. The current example has two regions which range from channel 5 to 163 and from channel 270 to 980. Currently the baseline regions cannot be defined by the cursor using the SSA.
- basefit(191,1,3) Set the baseline region in channels (e.g, [x1:x2][x3:x4]....) >> [5:163][270:980] >> T
When a spectral line is broadened by thermal and turbulent motions the expected shape is Gaussian. Therefore a standard procedure is to fit the spectral profile with a Gaussian curve.
The DISH GUI can be used by selecting Gaussfit under the operations menu. The number of Gaussian's must be specified along with initial values for the peak intensity, the full-width half-maximum, and the center channel of the peak. These can also be set using the cursor.
The SSA function gauss can be used to fit a profile with a Gaussian at the command line. The user function gaussfit allows the current SD record (the spectrum displayed in the plotter) to be used or any specified SD record. Unlike the previous example user functions here we are not selecting a spectrum from the working set but assume that the SD record has already been selected.
For example, in the last section we fit and subtracted a polynomial baseline to scan 191, subscan 1. The following example fits a single Gaussian to one feature in this spectrum. Currently from the SSA only one Gaussian can be fit at a time and the cursor cannot be used to define the initial guesses.
- gaussfit() Enter the initial peak, height, and center (e.g., 25 2 185) >> 25 2 185 T
This chapter has described some basic AIPS++ data reduction techniques for use with the GBT. The ultimate goal is to provide the user with enough information such that end-to-end data processing can be performed. This document is expected to evolve as more experience is gained with the GBT.