Original: Juan Uson, 8/99
Updated: J. Hibbard, 7/05

VLA HI AIPS Data Reduction Guide

This is a step-to-step instruction how to calibrate, map, and reduce a typical observation of a small galaxy in the HI line using the VLA in the D-configuration. Some familiarity with AIPS fundamentals is assumed. Instructions are not necessarily complete; it is strongly advised to carefully read the HELP file for most of the tasks mentioned here. See also the AIPS COOKBOOK, especially Chapter 4: Calibrating Interferometer Data (ps.gz) and Chapter 8: Spectral-Line Software (ps.gz). For a more streamlined guide, see Appendix B: A Ste-p-by-step guide to Spectral Line Data Analysis in AIPS (ps.gz).

The two files needed are called SUM_SPEC_0.gz and SUM_SPEC_LINE.gz. Using ftp, transfer these files to your own workstation. Size of the files is 7.4 and 43 MB, respectively, so you may want to do this during a time of the day with low net traffic. Once on your own disk, these files have to be 'unzipped':

%gunzip *.gz

After this, the sizes are 15 MB and 70 MB.

Set up a variable pointing to the directory where you put the relevant files to be read into aips.

%setenv MYFITS '/users/myname/fitsdirectory/'
%aips

pick a printer, then enter your aips number Next, reset aips parameters and load in the procedures you'll need for the default calibration:

   >RESTORE 0
   >RUN VLAPROCS

Make it so that the command line comes back after starting each aips task.

   >DOWAIT=0

see how much free disk space there is:

   >FREE

Chose the disk with the largest free space which isn't marked as RESERVED. If this is disk 'n' then designate it by:

   >INDISK n
   >OUTDIS n

Now read your data into aips:

FITLD or FILLM

FITLD is used to read in the UV data from the FITS files SUM_SPEC_0 and SUM_SPEC_LINE. FILLM is used for archival data or native For the FITLD case:

    >TASK 'FITLD'   
    >UVCOMP=0 
    >CLRONA 

This makes OUTNAME, OUTCLASS etc. blank

    > GO

This runs the task with the given inputs. In this case, The datasets will both have NAME='13/04/95'; one has CLASS 'CH0' and one class 'LINE'. In subsequent steps, the flux and phase calibration will be done on the CH0 data and then transferred to the LINE data.

For the FILLM case:

    >TASK 'FILLM'
    >CLRONA 
.   >UVCOMP=0 
    >DOUVCOM=-1
    >DOWEIGHT 2
    >CPARM 0;CPARM(7)=150
    >GO

    >UC

Gives you a listing of the UV data in your user area.

For VLA Archival data:

Example:
want to read in 2 files VLAOBS:
   AV220_1
   AV220_2

   >INFILE 'T:AV220_'
   (where T is you export defined directory that contains both files)

for the first file, AV220_1

    >NFILES 0
    >NCOUNT 1
   (this starts at 0 and increments to file _1 so to read in AV220_1)
for the second file:
    >NFILES 1
    >NCOUNT 1
  (to read in AV220_2  ..and so on )

    >GETN <ch0>;zap

We don't want to use the on-line "channel 0". We will make a new version below using AVSPC.

    >GETN <line UV data>
    >IMH

Look carefully at the header. Check the number of frequency channels (127). RA and Dec are 0 since there are several sources in the dataset (multi source file). A list of the sources is in the SU (source) table. You can look at any table using PRTAB. Note the number of visibilities in the CH0 data and the LINE data: they should be identical. In CH0, each visibility contains just one spectral line channel; in LINE, there are 127 channels. Also note the listing of tables attached to the dataset. See CH1 of the AIPS cookbook for a description of these tables.

In the following calibration steps, AIPS simply writes values into tables. It does not modify the raw UV data at all until the SPLIT step. To undo calibration or flagging, simply delete the relevant calibration (SN or CL;2 tables) or flagging (FG) tables as follows:

    >INEXT 'CL'
    >INVER 2
    >EXTD

BE VERY CAREFUL TO NOT EVER DELETE CL TABLE VERSION 1! IF BY ACCIDENT YOU DO, YOU WILL HAVE TO READ IN YOUR DATA FROM SCRATCH

The calibration procedure is:

• SETJY and GETJY write flux values into the SU table
• VLACALIB writes calibration coefficients into the SN table for times when the array was observing calibration sources
• VLACLCAL interpolates calibration coefficients from an existing SN table at regularly spaced time intervals, writing them into CL table version 2 (CL table version 1 contains 1's for all antenna gains and 0's for phases)
• SPLIT applies the CL table version 2 to the UV dataset

PRTAN

Get a printout of antennae positions, for use in flagging and chosing a reference antenna:

    >TASK 'PRTAN';DOCRT 0;GO

LISTR

Get a printout of the observing sequence, frequencies, etc.

    >TASK 'LISTR'; 
    >OPTYPE='SCAN'; DOCRT=-1; OUTFIL=''; GO

Pick up the output and look over the observing sequence. Identify the flux calibrator, phase calibrator, source, and observing sequence. Also the total time on each source. If the LISTR output shows that you have more than one FREQID in the observation, copy each out to its own file. If only one FREQID is present (as it is in the example dataset) you can skip the next step. In the practice dataset, the flux calibrator is 1331+30500002, the phase calibrator is 1445+09900002, and the source is N5921. There is only 1 freqid (corresponding too...) The flux calibrator was observed... for ... ** and the phase calibrator was observed... for ... ** The source was observed... **

UVCOP [may not be necessary]

Make a dataset containing only a single frequency

    >TASK 'UVCOP';
    >FREQID=n;OUTNA=<give name>;OUTCL 'UVLINE';GO
    >GETN <uvcop output>;GO INDXR
    >GO 'LISTR'

This copies the data out for FREQID=n, re-indexes the observe sequence to include only sources observed with that FREQID, and then re-runs listr to get a new scan list output. Use that from now on. Repeat for each FREQID.

SETJY

The flux of the flux calibrator is well tabulated and its frequency dependence is known. SETJY will calculate this flux based on the frequency it finds in the header, and puts the result in the SU table.

   >TASK 'SETJY'
   >SOURCE '<flux calibrator>' ' ';OPTYPE 'CALC';GO

You can check that it was put into the SU table via:

   >INEXT 'SU';DOCRT 132;GO PRTAB

AVSPC

Create a data file that averages data over all channels except those at the band edges. For calibrators, this file will have a higher signal to noise than the channel data. It will be used to derive calibration and flag table entries, which will then be copied to the line data. The range of channels should be as many as you expect to use in your final datacube. They are specified by ICHANSEL, where the syntax is ICHANSEL [start ch] [stop ch] [increment] 1. In practice, these should be as follows:

• 127 channel dataset: ICHANSEL 5 124 1 1
• 63 channel dataset: ICHANSEL 4 59 1 1
• 31 channel dataset: ICHANSEL 3 29 1 1

   >TASK 'AVSPC'
   >ICHANSEL 5 124 1 1
   >OUTNA INNA; OUTCL 'UVCH0';GO

TVFLG [calibrators, no calibration]

Run TVFLG on the CH0 data, only looking at the calibrator scans (SOUR='' and CALCO '*'). In this case, the data are of sufficient quality that flagging of bad data is not necessary, but TVFLG is an excellent way to get acquainted with the data.

   >TASK 'TVFLG'
   >GETN <ch0>
   >DOCAL 0; DOCAT 0;SOUR ' ';
   >CALCO '*'   
   >DPARM 0;DPARM(6) <integration time>;GO

You should know the intgegration time from the observe file. If it is archival data, you may have to figure it out via trial and error. If the integration time is set too short, blank lines will appear in between each line of real data. Change DPARM(6) until the data are displayed without blanks, except for the times between calibrator scans.

Most useful ways to dislay data here: delta amplitude; delta vector amplitude. Only flag grossly incorrect things: bad first scans; obviously bad times (horizontal features); obviously bad baselines (vertical features).

When flagging, make sure that the "all channels" flag is set (see writing at bottom of screen). This way, when you copy the FG table from your CH0 data to the LINE data, the flags will apply to your channels.

You also probably want to make sure the "all source" flag is set, so that when you flag an antenna or baseline it will apply also to your source data. But then you must be careful if you select the "flag timerange" option to not flag in the area between calibrator scans (unless you really mean to).

You will only want to flag the stokes data that you are looking at. To flag RR only, set stokes flag to 1000. To flag LL, set stokes flag to 0100. To flag both RR and LL, set stokes flag to 1100.

VLACALIB

Begin of the actual calibration. First, need to know if there are any UVLIMITS for either of your calibration sources. These are the range of spatial frequencies over which the calibrators look like unresolved ("point") sources. You get this information from the VLA calibrator manual, which can be found online via the following navigation: www.nrao.edu->VLA->Astronomers->Under "VLA Calibration": VLA Calibration Manual->Item 4.2 List of Calibrators (URL in 2005 is http://www.aoc.nrao.edu/~gtaylor/csource.html). Each source has two names, one based on B1950 coordinates, the other based on J2000. If the UVLIMITS do not overlap, you can run VLACALIB once with one set of limits and calibrate both sources. If they differ, you will need to run VLACALIB separately for each source.

In the practice dataset, the flux calibrator 1331+305 has the following entries in the calibration for L-band:

   BAND        A B C D    FLUX(Jy)    UVMIN(kL)  UVMAX(kL)
   =====================================================
    20cm    L  S S P P      15.00                       visplot

So it is a Prefered calibrator when the VLA is in D-array, with a flux at L-band of ~15 Jy, and no UV limits. The visplot link will plot its visibilities - print that out for your notes. Similarly, the phase calibrator 1445+099 is a Prefered with no uvlimits and a flux of ~2.6 Jy. In this case, there will be no UVLIMITs, and a single pass of VLACALIB with CALSOU ' ' and CALCO='*' will work.

    RUN VLAPROCS

   >TASK "VLACALIB'
   >CALSOU, UVRA, CALCO=<values based on calibration manual>
   >MINAMP 10; MINPHS 10; (only print out closure errors that
             deviate by more than 10% or 10deg)
   >REFANT n (chose an antenna near the center of the array, as
              determined from the PRTAN output)
   >DOCAL 0;FLAGVER 1;DOPRINT 1
   >GO VLACALIB

Three things will be printed out: closure errors, table of antenna gains for times observing each calibrator, and table of antenna phases (relative to the REFANT). There will be separate solutions for each polarization and each IF. The closure errors should not be too far from MINAMP, MINPHS. From the tables, check that the amplitudes are reasonably constant and that there are no large phase jumps. If they do, identify and flag the bad data using a combination of TVFLG, UVPLT, UVFND, UVFLG. It is often easier to find bad data after calibration has been applied, so at this step I usually only look for obviously bad baselines, times or antennas.

GETJY

Based on the flux of the flux calibrator and on the values in the two SN tables, GETJY will calculate the flux of the phase calibrator and put it in the SU table. This process is necessary since the fluxes of phase calibrators are not as constant and well monitored as are the few flux calibrators in the sky. Again, check with PRTAB on the SU table that the value found has been put in the table.

   >TASK 'GETJY'
   >SOUR '<phase cal>'
   >CALSO '<flux cal>'
   >GO

Keep note of the reported value.

VLACLCAL

Interpolates SN table(s); puts results in CL table version 2. A 'skeleton' CL table was already present in version 1. How you interpolate the values in the SN table depends on the observing sequence. If each time on the source is preceded and followed by a scan on the phase calibrator, then use SIMP interpolation. Otherwise, use a BOXcar. First, calibrate your source by interpolating the SN solutions from the phase calibrator:

   >TASK 'VLACLCAL'
   >REFANT n
   >SOUR '<source>'
   >CALSO '<phase cal>'
   >INTERPOL 'SIMP' (or INTERPOL 'box')
   >BPARM 0 (or BPARM 1,1,0 for INTERPOL 'box' and a 1hr boxcar)
   >GAINVER 0;GAINUSE 2;GO VLACLCAL

Next, calibrate your calibrators by using themselves:

   >SOUR ' ';CALSO ' ';CALCO '*';GAINUSE 0
   >GO

UVPLT

Check calibration on phase and flux calibrator

   >TASK 'UVPLT'
   >UVRA 0; source <phasecal>
   >XINC 5;
   >DOCAL 2;GAINUSE 2;
   >APARM 0; DOTV 1;GO

If this takes too long, then decrease XINC. If the plot looks too sparse, increase it. Make a hardcopy of the results for your records:

   >DOTV 0;GO;WAIT
   >DPARM 0;GO LWPLA

All visibilities should have uniform value of level specified by SETJY or GETJY, at least over applicable uvlimits. If not, something is seriously wrong. You will need to track down the bad data using tasks like TVFLG, UVPLT, UVFND, UVFLG, then delete the current calibration (SN table and CL table version 2: INEXT 'SN';INVER 0;EXTD;INEXT 'CL';INVER 2;EXTD) reset the flux values in the SU table (SOUR ' ';OPCO 'REJY'; GO SETJY), and redo all steps from SETJY to UVPLT.

TVFLG [calibrators, with calibration]

Rerun TVFLG on calibrators (CALCO '*'), this time with DOCAL=2. All baselines should have about the same amplitude now. Besides the displays you used last time, its also good to look at amplitude and delta phase.

IMAGR [phase cal ch0]

If UVPLT looks good, make sure there are no major imaging errors (see chapter on Error Recognition in Synthesis Imaging book) by mapping your phase and/or bandpass calibrator. IMAGR is an all-in-one task which applies the complex gains in the CL and BP tables (i.e. performs the actual calibration), does the Fourier transform to make maps in the sky plane, and optionally "cleans" the data as well. Map CH0 of the phase calibrator first; the result should be a point source in the field center, just like the beam. For cellsize, you want ~3 pixels across a beam width. The beam width is printed to the message window when you run IMAGR. Rough cellsizes: 12", 4", 1", 0.25" for data taken with the VLA in the D-, C-, B-, or A-arrays at L-band. Here, the image size doesn't need to be large, and should be a power of two. Try IMSIZ 64 64. Use ROBUST=-1 for maximum resolution.

   >TASK 'IMAGR'
   >OUTNA 'PHCAL.DR-1'
   >DOCALIB=2; FLAGV=1
   >ROBUST=-1
   >CELLSI 12 12
   >IMSIZE 64 64
   >NITER 500
   >DOTV 0
   >GO
   >MC

MC gives a listing of maps in your catalog. You should see the map you made (extension IIM), and the dirty beam (extension IBM). Get the map and display it to the tv:

   >GETN <phasecal map>
   >TVALL

FFT (includes COMB & UVNOU)

If the above has artefacts, you have to track them down. A very good diagnostic is looking at the baseline amplitudes and phases in the UV plane. You do this by doing a Fast Fourier Transform (FFT) of your map.

Running the task FFT on an image produces a Real & Imaginary FFT, which is basically all the antenna baselines as seen from the source. Short baselines are in the middle (note the hole in the center, since there are no baselines with zero spacing), and long baselines are in the outer regions.

To make an amplitude and phase FFT:

   >TASK 'FFT'
   >OPCO ' '
   >OUTNA INNA;GO
   >TASK 'COMB'
   >GETN <real part of FFT>
   >GET2N <imaginary part of FFT>
   >OPCOD 'POLI';GO
   >OPCOD 'POLA';GO
   >GETN <poli output>
   >TVALL

Bad baselines should show up very clearly. You need to find them and flag them. Use either TVFLG, UVPLT or UVFND. If you flag significant data, it may be necessary to redo your calibration. (See steps outlined above under UVPLT notes).

Look at R. Bracewell's "The Fourier Transform and its applications" to help you understand what you might expect in the Fourier plan for various types of sources (gaussians, point sources, etc). Or there is a great applet at http://www.narrabri.atnf.csiro.au/astronomy/vri.html. So, e.g. a point source at the phase center should the same amplitude on all baselines, and no phase structure (random phases). A point source offset from the phase center should have a uniform amplitude and a sinusoidal phase pattern. Try some at the above applet.

This is a really powerful diagnostic. The reason is that real sources, or weak sources, NEVER look like the array or are confined to a discrete number of points in the "UV" plane (which is basically what the FFT is showing you). So any high points that are confined to a few number of points is showing you problems with a correlator (if only 2 points - the UV point and its complex conjugate), or a baseline (if an arc), or a set of points (e.g. cross talk, which affects adjacent antennas). It also shows you when you have problems with short baselines or long baselines.

Sometimes its good to FFT the continuum, since the S/N is highest there, but often its best to FFT a continuum subtracted channel with little or no signal, since that should have NO structure in either amplitude or phase. But you don't have spectral line data, so I just throw this out there FYI. If you expect the source to be unpolarized, then you can use the stoke V, U maps as good diagnostics.

The trouble is tracking down the offending data in the UV visibilities based on the FFT map - there is no easy way to go between them in AIPS. Use TVFLG, UVFND, UVPLT. For example, if you see e.g. a strong SW-NE stripe in your map, this suggests a problem with a NW-SE baseline. After an FFT, you might see are a couple of whopping points, which should be one baseline at one time with a high amplitude. This should be fairly easy to find using TVFLG (see below), or even by doing a UVPLT using the defaults (plots amplitude vs. baseline length; set XINC equal to something like 2-5 and DOTV=1; if not getting enough points then decrease XINC) - the vast majority of points will cluster at some reasonable flux level, and there will be a few at high amplitude. These need to be clipped (either CLIP or TVFLG), and then re-image, re-FFT and see if there are still problems at this position angle. Another example is if you see an arc of high amplitudes in the FFT - this means a single baseline is bad for a certain period of time. If you see a set of arcs, it means an antenna is bad for a certain period of time.

A special case is when you find the high points all lying along the x-axis in the FFT. This is a sign of RFI (see the EXPLAIN for the task UVNOU). To get rid of it, use the task UVNOU. The format for UVRANGE is different for this task than most others; instead of it being an annulus in UVdistance, its the distance away from the U and/or the V planes. You want to use UVRANGE= [a small number],0. I determine the small number by trial and error - start with something like 100 and increase or decrease based on what you see in the FFT. You want to flag just the points very close to the axis. NOTE: unlike some other flagging tasks, UVNOU flags the actually UVdata instead of producing a flag table. So if you run UVNOU on your ch0 data, you need to run it with the same parameters on your line data.

When satisfied:

TVFLG [source, with calibration]

Now run TVFLG on your calibrated source data

   >TGET TVFLG
   >CALCO ' ';SOURCE 'N5921';DOCAL=2;GAINUS 2;FLAGVER=1;GO

In general, this will not look like your calibrator data (similar amplitudes on all baselines). Set "Sort by length" - usually, your shortest baselines (now displayed to the left side of the tv) will have more flux than your longer baselines. Flag stuff that sticks out in time or baseline. After that, "Enter Scan Time" of about 10 integration times (performs a rolling boxcar average), and display amplitude then phase differences, and clip large deviations ("clip interactively). Note that no matter how much you flag here (and it shouldnt be a lot), it does not affect your calibration, so you will not need to delete tables and redo anything.

IMAGR [source ch0]

Now map the continuum sources in the field of your source. Map the entire primary beam (30' FWHM at L-band). Set DOTV 1 to run TVFLG in interactive mode. Set boxes tightly around any sources, update them as cleaning progresses.

   >TASK 'IMAGR'
   >OUTNA 'N5921D.CH0-1'
   >DOCALIB=2;FLAGV=1;ROBUST=-1
   >ROBUST=-1
   >CELLSI 12 12
   >IMSIZE 512 512
   >NITER 500
   >DOTV 1
   >GO

If problems, try FFT/COMB, UVPLT, etc.

If ok, then spatial calibration is done; copy CL;2 and FG tables to the line data, and proceed with Bandpass calibration.

TACOP

Copy table CL version 2 and FG table from CH0 to LINE data.

   >TASK 'TACOP'
   >GETN <line uvdata>
   >outna inna;outcl incl;
   >GETN <ch0 uvdata>
   >INEXT 'fg';inver 1;outver 1;go
   >INEXT 'CL';inver 2;outver 2;go

This assumes that the complex gains can be split in a time dependent part and a frequency dependent part. CALIB and CLCAL take care of the time dependent - but frequency independent - part; in the next step we will determine the frequency dependent complex gains.

BPASS

Creates BP table; gains across band. Use 1331+30500002 as bandpass calibrator (it is the only calibrator which is bright enough that it has sufficient S/N in each channel to do a proper bandpass calibration). There are many ways to do this. The following calculates the bandpass response by dividing each visibility by its average over the many channels (as specified by ICHANSEL).

   >TASK 'BPASS'
   >GETN <line uvdata>
   >CALSO '<bpass calibrator>'
   >UVRA 0; DOCAL 0;FLAGVER 1;SOLINT 0;DOBAN 0;
   >REFANT n (same as vlacalib)
   >BPASSPRM 0,1,0,0,-2,0
   >ICHANSEL <same as AVSPC>
   >GO

POSSM

Plot the bandpass response for each antenna, and make sure no sudden jumps, which would indicate RFI.

   >TASK 'POSSM' 
   >APARM 0,1,0.9,1.1,-10,10,0,2,0
   >NPLOTS 9
   >DOTV 1

Inspect the bp table for each antenna, showing 9 plots per page. If ok, make hardcopies for posterity.

   >DOTV -1;go

Now check how flat your bandpass solution makes your phase calibrator calculating the global spectrum of your phase calibrator with both bandpass and calibration applied

   >APARM 1,0
   >NPLOTS 0
   >SOUR '<phase cal>'
   >DOCAL 2; DOBAND 1
   >DOTV 1

Should be flat, and at the flux level given by GETJY. If ok, make hardcopies for posterity.

   >DOTV -1;go

Print out copies of the above plot files. Do an >IMH on your uvdata to find the number of plot files (call this n), then:

   >task 'lwpla';dparm 0
   >for i=1:n;plver i;wait;go;end 

If you only observed your BP calibrator once, you are done. If you observed it several times, you have another decision to make - how to apply the bandpass. Type ">help doband" for the options. To help decide, repeat the last POSSM option with DOBAND=2, and again with DOBAND=3. Print out each. Pick the option which makes the phasecal spectrum look flattest. Quantify the flatness of the bandpass by looking at the maximum deviation across the band. Call the devaition delta over N channels, then the bandpass is flat to Delta/Flux/Nchan. This is a good quantity to know.

ALTDEF

For older data, it is necessary to enter the velocity which corresponds to the observed sky frequency relationship into the header. First do an IMH on the UVData. If it does not have lines that look like the following:

AIPS 1: Rest freq   1420.406         Vel type: OPTICAL wrt SUN
AIPS 1: Alt ref. value  1.4800E+06  wrt pixel   32.00

then you will need to do this. You need to know the observed velocity either from the original LISTR 'scan' output or the observe file. For the practice data:

    >AXTYPE 'OPTHEL'
    >AXVAL=<the radial velocity of source in m/s>
    >AXREF=<central channel>
    >RESTFREQ=1420E9,40575
    >ALTDEF

IMAGR [phasecal channels]

It is good practice to make a channel cube of the phase cal first; the result should be a point source in the field center, just like the beam. I always use BCHAN=1, so I can easily go from a plan in the image cube to the correspinding channel in the dataset. ECHAN=125 to avoid the band edges. Clean the image lightly, i.e. NITER=200.

   >TASK 'IMAGR' 
   >OUTNA 'phcal<array>-1'
   >DOCALIB=2; GAINUSE 2; FLAGV 1;
   >DOBAND=n (from POSSM) 
   >BCHAN 1; ECHAN 125
   >CELLSI 12 12
   >IMSIZE 64 64;
   >UVWT=''; ROBUST=-1; 
   >NITER 200
   >DOTV 0

View using TVMOVIE, LTYPE=6 so that the channel number is displayed in each frame. You might have to change PIXR in order to see the signal, e.g.

   >GETN <phasecal cube>
   >PIXR -0.01,0.03
   >TVMOV

ISPEC

Check that spectra are flat using ISPEC:

   >TASK 'ISPEC'
   >GETN <phasecal data cube>
   >OPTY 'FLUX';LTYPE 6
   >BLC 0 0 4;TRC 0 0 125;DOTV 1;GO 

This better be flat with the flux reported by GETJY Make a hardcopy (dotv -1;go;plver 0;go lwpla).

IMAGR [source line+continuum data]

If that looks good, re-run IMAGR on your source with the continuum still in it. We will not keep this; its to get an idea of which channels contain line emission. Note that IMAGR has adverbs BCHAN and ECHAN, so it needs to be run only once to do all channels. You may want to first run IMAGR on a single test channel (e.g BCHAN=64; ECHAN=64); this also allows you to test your mapping and cleaning parameters. Use a IMSIZE large enough to map the primary beam (30'). CELLSIZE as before. Shallow clean to remove large sidelobes.

   >TGET IMAGR
   >SOURCE='N5921'
   >OUTNA 'N5921D.r=1c'
   >ROBUST=1
   >IMSIZ 512 512
   >NITER 500;DOTV 0;
   >GO

The result will be a 'cube' of maps and a 'cube'of beams, with frequency along the third axis. Check this with IMHEAD. Make note of the "Conv size" & "Position angle" from the header - this is the size (FWHM) of the Gaussian used to replace dirty point sources in the original map. It is therefore the resolution of your map. You will need it below. This map contains signal from both continuum and spectra line sources. If it looks good, then its finally time to split the source data into its own dataset, with all calibration and flags from the tables applied to the data. TVMOV with LTYPE=6 a few times and try to identify the first and last channels where line emission appear.

SPLIT

Use SPLIT to convert your multisource file into datafiles containing just a single source. SPLIT applies the complex gains in the CL and BP tables (i.e. it performs the actual calibration).

   >TASK 'SPLIT'
   >SOUR 'N5921'
   >CALCO ' ';DOCAL 2;GAINUS 2;FLAGVE 1
   >APARM 0;NCHAV 1;ICHANSEl 0
   >DOBAND=<from POSSM>;
   >GO

I'll generally rerun the last IMAGR again, using the newly split dataset (but now setting DOCAL=-2; DOBAND=-1), to make sure everything got applied correctly.

Give the SPLIT file a descriptive name:

   >GETN <split dataset>
   >OUTNA 'N5291D';OUTCL 'UVLINE';RENA

---

OPTIONAL: Self Calibration

• AVSPC line free channels from split data
• IMAGR Robust-1, UVRANG, DOTV 1, use tight boxes
• CCMRG on resulting image
• PRTCC to find number of clean components to use
• CALIB getna [AVSPC data], get2n [AVSPC map], ncomp=from prtcc, uvra as in IMAGR, calso [source], solmo 'p' for phase only selfcal, refant [as in original vlacalib], solint 10, aparm 4,0,0,0,2,0, cparm 0, go
• IMAGR new uvdataset. If better, EXTD 'sn' file from avspc data, ccmrg,prtcc,calib, etc.
• TACOP SN file from AVSPC data to SPLIT UVLINE data
• SPLIT UVLINE data with SN table applied. OUTNA .SC

---

TASAV

Make a dummy file which contains no data, but all the calibration tables. You will keep the split data and resulting products, but there is no need to keep the raw data, since you can always get this back from the archive. But if you keep a file with your calibration tables attached, then if you do need to go back to the raw data, you can easily copy the calibration and flags from this file.

    > Task 'TASAV'
    > GETN <uvdata before split>
    > OUTNA 'N5921D.MS';OUTCL 'TASAV';GO

ZAP

This is a good time to get rid of files you no longer need. Basically, at this stage you only really need the SPLIT and TASAV uvdata, and the best cleaned ch0 map. We will make a new ch0 data later. We don't need any phase calibrator maps, nor source maps with continuum in them. So for each file you don't want:

    >GETN <i>
    >ZAP

To do a whole sequence of files, set up a loop:

    >FOR i=4:18;getn i;zap;end

When this is cleaned up, get remaining data into adjacent slots on disk: >recat

UVLIN

Now its time to remove the continuum from the line data. Subtracting a continuum from UV data can be done using the task UVLIN. This is an iterative procedure: UVLIN it, map it, review it; if negatives appear or ISPEC doesn't look right, try a new set of channels and repeat.

From the shallowly cleaned source datacube, you identified line-free channels. Lets say that the line emission appears from channel A to B. We'll guess that the line free channels are b=A-2, c=B+2. For the test data, these are ???? We also need the first and last channels used in our original AVSPC (used to make ch0). Call these a and d (e.g., for the test dataset these were 5, 124. So the line-free channels that we want to fit are a->b and c->d. These are specified in UVLIN using ICHANSEL=a b 1 0 c d 1 0

But before running, we note that UVLIN also flags data that deviate in some channel by more than FLUX from the mean over the included channels. See EXPLAIN UVLIN for more on this. In practice, I determin FLUX interactively - I start with FLUX~3, and increase or decrease it based on how many correlators it flags - on otherwise well behaved data, I want to flag ~<1% of the correlators.

   >TASK 'UVLIN' 
   >GETN <UVLINE DATA>
   >OUTNA INNA
   >OUTCL 'UVLIN'
   >DOCONT 0;ORDER 1;PRTLEV 1
   >ICHANSEL=5 ?? 1 0 ?? 125 1 0
   >FLUX 3
   >GO

IMAGR,ISPEC,UVLIN, repeat

Make shallowly cleaned datacube, inspect it, use ISPEC to take the spectra through the location of continuum sources (get the pixel location of continuum sources from your ch0 map made with the same CELLSI, IMSIZE as your channel map using the verbs CURV or IMPOS, and use these for BLC,TRC in ISPEC) and your line source. When you are satisfied that UVLIN only included line-free data, make a deeper cleaned datacube.

IMAGR [UVLIN'd source, R=+1]

Use the same inputs as before, but a larger NITER Also set BMAJ,BMIN,BPA so that each channel has the exact same resolution. Use Robust=1 to improve sensitivity to extended resolution. ??(NBOX? CLBOX?)

   >TGET IMAGR
   >SOURCE='N5921'
   >OUTNA 'N5921D.r=1c'
   >ROBUST=1. 
   >BMAJ <major axis FWHM from previous IMAGR run, rounded off>
   >BMIN <minor axis FWHM from previous IMAGR run, rounded off>
   >BPA <position angle from previous IMAGR run>
   >IMSIZ 512 512
   >NITER 1000;DOTV 0;
   >GO

Get the cube and change the 3rd axis from frequency to velocity, using ALTSW Also, I always rename cleaned cubes with the OUTCL 'LMV', which designates the cube order as being RA, DEC, VEL.

   >GETN <r+1 cube>
   >ALTSW
   >IMH (notice units are now m/s on 3rd axis, and increment is
         your channel spacing in m/s)
   >OUTNA INNA,OUTCL 'LMV';RENAM

IMAGR [UVLIN'd source, R=-1]

I always also make a R=-1 version of the data. You'll have to run it on the central channel first to get BMAJ,BMIN,BPA

   >TGET IMAGR
   >ROBUST=-1
   >BCHAN 64
   >ECHAN 64
   >OUTNA 'ch64.r-1'
   >GO
   >GETN <ch64.r-1>
   >IMH (make note of beam parameters)
   >TGET IMAGR
   >BCHAN 1
   >ECHAN 125
   >OUTNA 'N5921D.r=1c'
   >GO

ALTSW and rename as above

   >GETN <r+1 cube>
   >ALTSW
   >OUTNA INNA,OUTCL 'LMV';RENAM

When satisfied, we want to make two new ch0-like datasets from the split (i.e. not UVLIN'd) datacube. _1 A new ch0, but this one should be nicer b/c the bandpass shape has been removed and we'll use UVLIN to flag outliers _2 A line-free continuum map, made from line free channels.

UVLIN (flagging only)

First run UVLINE only to flag outliers; retain the continuum

   >TGET UVLIN
   >DOCONT 1
   >OUTCLA 'LINE+C';GO

AVSPC

now make the new ch0 and line-free continuum

   >TASK 'AVSPC'
   >GETN <uvlin docont=1 output>
   >ICHANSEL 5 125
   >OUTNA INNA
   >OUTCLA 'UV CH0'
   >GO

   >INCHANSEL 5 *A* *B* 125
   >OUTCL 'UVCONT'
   >GO

IMAGR [new Ch0, line-free continuum]

Use IMAGR to map both of those datasets. Use same CELLSI, IMSIZE so that you can correlate continuum sources with line data.

   >TGET IMAGR
   >ROBUST -1
   >GETN <ch0>
   >OUTNA 'N5921Dr-1.ch0'
   >NITER <large number - depends on number of sources. ~5000>
   >DOTV 1
   >GO

   >GETN <uvcont>
   >OUTNA 'N5921Dr-1.nol'
   >GO

IMEAN

Measure noise in each map using TVWIN, GO IMEAN Make note of this value. Also note number of pixels per beam. These are needed later.

TRANS

Many routines that work on cubes and velocity profiles require the data with frequency - or velocity - on the first axis. This is done using the task TRANS with TRANSCOD='312' on the cleaned cube. Use IMHEAD on the result to check whether TRANS did what you wanted.

   >TASK 'TRANS'
   >GETN <r+1 datacube>
   >BLC 0; TRC 0
   >TRANSCO '312'
   >OUTNA INNA
   >OUTCL 'VLM'
   >GO

   >GETN <r-1 datacube>
   >OUTNA INNA
   >GO

---

JEH...I'M HERE...

MOMNT

MOMNT creates maps containing the various moments of the velocity profiles. It works on the cube you just created with TRANS. Read the HELP file for MOMNT carefully, and play around with the values of CELLSIZE, ICUT, and FLUX.

   >TASK 'MOMNT'
   >CELLSI 3 5
   >BLC a 0 0
   >TRC b 0 0
   >OUTCL '0'
   >FLUX -1
   >ICUT (1sigma)
   >GO


 # ICUT~ 1-1.5* σ
 #repeat for ZMOM0: 

            >FLUX -1
            >ICUT~ 1-1.5* σ (Jy/b)
            >OUTN INNA
            >OUTCL '0'
            >GO
 
            MOM0,1,2 
            >FLUX 0
            >OUTN INNA
            >OUTCL '012'
            >GO  

-- AmandaH - 25 Jul 2005

---

NOTES:

Naming convention:

I find it very useful to give each file a descriptive name so that you can tell what it is without having to do an IMH. My naming convention is:

UVDATA:
OUTNA='<source><array><MS if dataset contains multisource data>'
OUTCL=one of the following:
   'UVLINE': uv line data.
   'UVCH0': average over many channels.
   'UVCONT': average over line-free channels.
   'UVLIN': uv line data after UVLIN.
   'TASAV': uv dummy file with all calibration tables attached.

MAPS:
OUTNA='<source><array>.<Robust value><ch0 or cont or nothing>'
OUTCL=one of the following: 
   'ICLN': cleaned single channel map.
   'IMAP': uncleaned map or cube.
   'LMV': cleaned cube in [x,y,z]=RA,DEC,Vel order
   'VLM': cleaned cube in [x,y,z]=Vel,RA,DEC order
   'ZMOM0': mom0 made with flux=-1
   'MOM0,MOM1,MOM2': moment maps made with flux=0

To change the name of a file:

>getn <file no.>
>outna 'n5921d.ms'
>outcl 'uvline'
>renam

---

Final Products:

Print Outs:

• Antenna positions from PRTAN
• Original LISTR OPTY='SCAN'
• Final LISTR OPTY='SCAN' (after UVCOP)
• printout from original VLACALIB
• printout from final VLACALIB

Plots:

• UVPLT for each of your calibrators and source
• BP table for each antenna (POSSM, APARM(8)=2)
• Spectrum of Phase Cal with BP applied (POSSM, APARM(1)=1)
• ISPEC of source
• ISPEC through cube at location of brightest continuum source
• KNTR of line-free continuum greys with Ch0 contours
• KNTR of channel maps (contours on optical greys if possible)
• KNTR of mom1 contours on greys of mom0
• KNTR of mom1 contours on greys of mom2

UV data:

• [name][array].MS.TASAV (dummy uvfile with BP,CL,SN,FG tables)
• [name][array].UVLINE (uv line data)
• [name][array].UVCH0 (average over many channels)
• [name][array].UVCONT (average over line-free channels)
• [name][array].UVLIN (uv line data after UVLIN)

MAPS:

• [name][array]-1.CH0.ICLN (final ch0)
• [name][array]-1.-L.ICLN (line-free continuum)
• [name][array].R-1.LMV
• [name][array].R-1.ZMOM0
• [name][array].R-1.MOM0
• [name][array].R-1.MOM1
• [name][array].R-1.MOM2
• [name][array].R+1.LMV
• [name][array].R+1.ZMOM0
• [name][array].R+1.MOM0
• [name][array].R+1.MOM1
• [name][array].R+1.MOM2

Things to document for each map:

• theoretical and measured noise (mJy/b)
• clean beam size FWHMxFWHM at P.A.
• For continuum maps: channels which went into making them
• For moment maps: icut, flux, cellsi,
• for moment 0 maps: total flux (Jy km/s)

---

Calculating Theoretical Noise:

 S_rms  =  k*fw/sqrt(N_ANT*(N_ANT-1)*N_IFs*Δt_hr*Δν_MHz)  [mJy/beam]

 k~8.0 @ L-band (from VLA "observational status summary")
 fw ~ 1.03 (Robust = +1)
    ~ 1.2-1.4 (Robust =-1 (see task 'imagr' output)
 N_ANT ~27
 N_IFs = # of polirization pairs (LL & RR)
 Δt_hr = time on source from observing log in hours
 Δν_MHz = bandwidth (from imh of data cube)

-- AmandaHeiderman - 25 Jul 2005

Measuring Noise:

>GETN <final image cube R=+1>
>TVWIN (select area near source)
>GO IMEAN

write down rms, # pixels/beam, and # pixels in chosen area
repeat for R=-1 cube

-- AmandaHeiderman - 25 Jul 2005

Measuring Flux:

Flagging based on UVSUB residuals

-- JohnHibbard - 07 Sep 2005
-- AmandaHeiderman - 08 Sep 2005 (updated)