From akemball@aoc.nrao.edu Sat May  5 20:11:22 2001
Date: Fri, 4 May 2001 17:30:40 -0600 (MDT)
From: Athol Kemball <akemball@aoc.nrao.edu>
To: aips2-naug@zia.aoc.NRAO.EDU
Cc: akemball@zia.aoc.NRAO.EDU
Subject: Spectral line script



I have appended a Glish script which reduces the H I test dataset on
NGC 5921, which is a designated test dataset on Steve's user group
page.  It was also used in the summer school last year. This script
does channel zero and bandpass calibration, with image plane continuum
subtraction (~ok in this case). I might suggest that someone refine
this next week, once our build transition is complete. A small fix for
the fact that these data only have RR,LL was required; this isn't in
stable, which is about 5 weeks old at the moment. We expect a new
stable next week once the build re-organization is complete.  Until
then daily and weekly remain unavailable as per my earlier message
this week.

We need a few more line datasets. The ones I have previously reduced
are not public, and can't be added to our test database. They also
need to be of modest size (few 10's MB so they don't fill up our
database area).

The resulting image can be found at:
      http://www.nrao.edu/~akemball/ngc5921_movie.gif

George converted the cube to the GIF animation.

-Athol

--------------------------------------------------------------------------

include 'imager.g';
include 'calibrater.g';
include 'image.g';
include 'ms.g';
include 'flagger.g';

#
# VLA H_I line data reducion: NGC 5921
#
# Load the data
#
   fitsfilename:='/export/home/voyager/akemball/fits/vla_line/LINE.UVFITS';
   m:=fitstoms(fitsfile=fitsfilename, msfile='ngc5921.ms');
   m.summary(verbose=T);
   m.done();
#
# Basic editing
#
   fg:=flagger('ngc5921.ms');
   fg.setflagmode('unflag');
   fg.query(query='FIELD_ID >= 0');
#
# Flag data for 1331+305
#
   fg.setflagmode('flag');
   fg.setids(fieldid=1);
   fg.filter(column='DATA', operation='range', comparison='Amplitude',
             range=['0.0Jy','1.6Jy']);
#
# Flag data for 1445+099
#   
   fg.setflagmode('flag');
   fg.setids(fieldid=2);
   fg.filter(column='DATA', operation='range', comparison='Amplitude',
             range=['0.0Jy','0.4Jy']);
#
# Flag data for NGC 5921
#
   fg.setflagmode('flag');
   fg.setids(fieldid=3);
   fg.filter(column='DATA', operation='range', comparison='Amplitude',
             range=['0.0Jy','0.11Jy']);
   fg.done();
#
# Set source model
#
   imagr:=imager('ngc5921.ms');
   imagr.setjy(fieldid=1, fluxdensity=[14.8009,0,0,0]);
   imagr.done();
#
# Solve for G Jones using pseudo-continuum
#
   cal:=calibrater('ngc5921.ms');
   cal.setdata(msselect='FIELD_ID <= 2');
   cal.setsolve(type='G', t=300.0, preavg=300.0, refant=15, 
                table='ngc5921.gcal');
   cal.solve();
#
# Establish the flux density scale
#
   cal.fluxscale(tablein='ngc5921.gcal', tableout='ngc5921.gcalout',
                 reference='1331+30500002',
                 transfer=['1445+09900002']);

   cal.plotcal(tablename='ngc5921.gcalout');
   cal.done();
#
# Solve for B Jones (bandpass)
#
   cal:=calibrater('ngc5921.ms');
   cal.setdata(msselect='FIELD_ID==1');
   cal.setapply(type='G', t=0.0, table='ngc5921.gcalout');
   cal.setsolve(type='B', t=86400.0, preavg= 86400.0, refant=15, 
                table='ngc5921.bcal');
   cal.solve();
   cal.plotcal(tablename='ngc5921.bcal');
   cal.done();
#
# Correct the line data for G and B Jones
#
   cal:=calibrater('ngc5921.ms');
   cal.setapply(type='G', t=0.0, table='ngc5921.gcalout', 
                select='FIELD_ID==2');
   cal.setapply(type='B', t=0.0, table='ngc5921.bcal');
   cal.correct();
   cal.done();
#
# Make a map of NGC 5921
#
   imagr:= imager('ngc5921.ms');
   imagr.setdata(fieldid=3, mode='channel', nchan=63, start=1, step=1,
                 spwid=1);
   imagr.setimage(nx=256, ny=256, cellx='10arcsec', celly='10arcsec', 
                  stokes='I', fieldid=3, spwid=1, start=1, step=1, nchan=63,
                  mode='channel');
   imagr.weight();
   imagr.make('ngc5921.model');
   imagr.clean(algorithm='clark', niter=3000, gain=0.1, threshold='0.0Jy',
               model='ngc5921.model', image='ngc5921.restored',
               residual='ngc5921.residual');
   imagr.done();
#
# Final restore
#
   imagr:= imager('ngc5921.ms');
   imagr.setdata(fieldid=3, mode='channel', nchan=63, start=1, step=1,
                 spwid=1);
   imagr.setimage(nx=256, ny=256, cellx='10arcsec', celly='10arcsec', 
                  stokes='I', fieldid=3, spwid=1, start=1, step=1, nchan=63,
                  mode='channel');
   imagr.weight();
   imagr.setbeam(bmaj='44.6arcsec', bmin='39.0arcsec', bpa='-174.1deg');
   imagr.restore(model='ngc5921.model', image='ngc5921.restored2');
   imagr.done();
#
#
# Subtract continuum channels (5-8) and (55-58) using image tool
#