##########################################################################  
#                                                                        #  
# Use Case Script for NGC 5921                                           #  
#                                                                        #  
# Converted by STM 2007-05-26                                            #  
# Updated      STM 2007-06-15 (Alpha Patch 1)                            #  
# Updated      STM 2007-09-05 (Alpha Patch 2+)                           #  
# Updated      STM 2007-09-18 (Alpha Patch 2+)                           #  
# Updated      STM 2007-09-18 (Pre-Beta) add immoments                   #  
# Updated      STM 2007-10-04 (Beta) update                              #  
# Updated      STM 2007-10-10 (Beta) add export                          #  
# Updated      STM 2007-11-08 (Beta Patch 0.5) add RRusk stuff           #  
# Updated      STM 2008-03-25 (Beta Patch 1.0)                           #  
# Updated      STM 2008-05-23 (Beta Patch 2.0) new tasking/clean/cal     #  
# Updated      STM 2008-06-11 (Beta Patch 2.0)                           #  
# Updated      STM 2008-06-13 (Beta Patch 2.0) demo version              #  
# Updated      STM 2008-06-14 (Beta Patch 2.0) post-school update        #  
# Updated      STM 2008-07-06 (Beta Patch 2.0) regression version        #  
# Updated      STM 2009-05-26 (Beta Patch 4.0) McMaster demo             #  
# Updated      STM 2009-12-02 (Release 0)                                #  
# Revised      JO  2010-04-13 (Release 3.0.1) edit viewer call           #  
# Checked      JO  2010-10-07 (Release 3.1.0)                            #  
# Updated      JO  2012-05-10 (Release 3.4.0)                            #  
#                                                                        #  
# Features Tested:                                                       #  
#    The script illustrates end-to-end processing with CASA              #  
#    as depicted in the following flow-chart.                            #  
#                                                                        #  
#    Filenames will have the <prefix> = ’ngc5921.demo’                   #  
#                                                                        #  
#    Input Data           Process          Output Data                   #  
#                                                                        #  
#   NGC5921.fits --> importuvfits  -->  <prefix>.ms   +                  #  
#   (1.4GHz,               |            <prefix>.ms.flagversions         #  
#    63 sp chan,           v                                             #  
#    D-array)           listobs    -->  casapy.log                       #  
#                          |                                             #  
#                          v                                             #  
#                     flagautocorr                                       #  
#                          |                                             #  
#                          v                                             #  
#                        setjy                                           #  
#                          |                                             #  
#                          v                                             #  
#                       bandpass   -->  <prefix>.bcal                    #  
#                          |                                             #  
#                          v                                             #  
#                       gaincal    -->  <prefix>.gcal                    #  
#                          |                                             #  
#                          v                                             #  
#                      fluxscale   -->  <prefix>.fluxscale               #  
#                          |                                             #  
#                          v                                             #  
#                      applycal    -->  <prefix>.ms                      #  
#                          |                                             #  
#                          v                                             #  
#                        split     -->  <prefix>.cal.split.ms            #  
#                          |                                             #  
#                          v                                             #  
#                        split     -->  <prefix>.src.split.ms            #  
#                          |                                             #  
#                          v                                             #  
#                    exportuvfits  -->  <prefix>.split.uvfits            #  
#                          |                                             #  
#                          v                                             #  
#                      uvcontsub   -->  <prefix>.ms.cont +               #  
#                          |            <prefix>.ms.contsub              #  
#                          v                                             #  
#                        clean     -->  <prefix>.clean.image +           #  
#                          |            <prefix>.clean.model +           #  
#                          |            <prefix>.clean.residual          #  
#                          v                                             #  
#                     exportfits   -->  <prefix>.clean.fits              #  
#                          |                                             #  
#                          v                                             #  
#                       imhead     -->  casapy.log                       #  
#                          |                                             #  
#                          v                                             #  
#                       imstat     -->  xstat (parameter)                #  
#                          |                                             #  
#                          v                                             #  
#                      immoments   -->  <prefix>.moments.integrated +    #  
#                          |            <prefix>.moments.weighted_coord  #  
#                          v                                             #  
##########################################################################  
print ’Demo Script for NGC5921 VLA HI observation’  
print ’Version for Release 0 (3.4.0) May-2012’  
print ’’  
 
import time  
import os  
 
scriptmode = True  
#  
# The prefix to use for all output files  
prefix=’ngc5921.demo’  
 
# Set up some useful variables (these will also be set later on)  
msfile = prefix + ’.ms’  
btable = prefix + ’.bcal’  
gtable = prefix + ’.gcal’  
ftable = prefix + ’.fluxscale’  
splitms = prefix + ’.src.split.ms’  
imname = prefix + ’.cleanimg’  
 
#  
# Get to path to the CASA home and stip off the name  
pathname=os.environ.get(’CASAPATH’).split()[0]  
 
# This is where the NGC5921 UVFITS data will be  
fitsdata=pathname+’/data/demo/NGC5921.fits’  
#  
# Or uncomment the following to use data in current directory  
#fitsdata=’NGC5921.fits’  
 
# Clean up old files  
# Use rmtables on ms and cal tables to clear cache  
# (not working on multiple runs for 2.4.0 release)  
#rmtables(msfile)  
#rmtables(btable)  
#rmtables(gtable)  
#rmtables(ftable)  
#rmtables(ftable)  
#rmtables(splitms+’*’)  
#rmtables(imname+’.*’)  
#rmtables(prefix+’.moments*’)  
 
# Final clean up of auxiliary files  
os.system(’rm -rf ’+prefix+’*’)  
 
#  
#=====================================================================  
#  
# Import the data from FITS to MS  
#  
print ’--Import--’  
 
# Safest to start from task defaults  
default(’importuvfits’)  
 
# Set up the MS filename and save as new global variable  
msfile = prefix + ’.ms’  
 
# Use task importuvfits  
fitsfile = fitsdata  
vis = msfile  
 
saveinputs(’importuvfits’,prefix+’.importuvfits.saved’)  
 
importuvfits()  
 
#  
# Note that there will be a ngc5921.demo.ms.flagversions  
# there containing the initial flags as backup for the main ms  
# flags.  
#  
#=====================================================================  
#  
# List a summary of the MS  
#  
print ’--Listobs--’  
 
# Don’t default this one and make use of the previous setting of  
# vis.  Remember, the variables are GLOBAL!  
 
# You may wish to see more detailed information, like the scans.  
# In this case use the verbose = True option  
verbose = True  
 
listobs()  
 
# You should get in your logger window and in the casapy.log file  
# something like:  
#  
# MeasurementSet Name:  /home/sandrock2/smyers/Testing2/Sep07/ngc5921.demo.ms  
# MS Version 2  
#  
# Observer: TEST     Project:  
# Observation: VLA  
#  
# Data records: 22653       Total integration time = 5280 seconds  
#    Observed from   09:19:00   to   10:47:00  
#  
#    ObservationID = 0         ArrayID = 0  
#   Date        Timerange                Scan  FldId FieldName      SpwIds  
#   13-Apr-1995/09:19:00.0 - 09:24:30.0     1      0 1331+30500002_0  [0]  
#               09:27:30.0 - 09:29:30.0     2      1 1445+09900002_0  [0]  
#               09:33:00.0 - 09:48:00.0     3      2 N5921_2        [0]  
#               09:50:30.0 - 09:51:00.0     4      1 1445+09900002_0  [0]  
#               10:22:00.0 - 10:23:00.0     5      1 1445+09900002_0  [0]  
#               10:26:00.0 - 10:43:00.0     6      2 N5921_2        [0]  
#               10:45:30.0 - 10:47:00.0     7      1 1445+09900002_0  [0]  
#  
# Fields: 3  
#   ID   Code Name          Right Ascension  Declination   Epoch  
#   0    C    1331+30500002_013:31:08.29      +30.30.32.96  J2000  
#   1    A    1445+09900002_014:45:16.47      +09.58.36.07  J2000  
#   2         N5921_2       15:22:00.00      +05.04.00.00  J2000  
#  
# Spectral Windows:  (1 unique spectral windows and 1 unique polarization setups)  
#   SpwID  #Chans Frame Ch1(MHz)    Resoln(kHz) TotBW(kHz)  Ref(MHz)    Corrs  
#   0          63 LSRK  1412.68608  24.4140625  1550.19688  1413.44902  RR  LL  
#  
# Feeds: 28: printing first row only  
#   Antenna   Spectral Window     # Receptors    Polarizations  
#   1         -1                  2              [         R, L]  
#  
# Antennas: 27:  
#   ID   Name  Station   Diam.    Long.         Lat.  
#   0    1     VLA:N7    25.0 m   -107.37.07.2  +33.54.12.9  
#   1    2     VLA:W1    25.0 m   -107.37.05.9  +33.54.00.5  
#   2    3     VLA:W2    25.0 m   -107.37.07.4  +33.54.00.9  
#   3    4     VLA:E1    25.0 m   -107.37.05.7  +33.53.59.2  
#   4    5     VLA:E3    25.0 m   -107.37.02.8  +33.54.00.5  
#   5    6     VLA:E9    25.0 m   -107.36.45.1  +33.53.53.6  
#   6    7     VLA:E6    25.0 m   -107.36.55.6  +33.53.57.7  
#   7    8     VLA:W8    25.0 m   -107.37.21.6  +33.53.53.0  
#   8    9     VLA:N5    25.0 m   -107.37.06.7  +33.54.08.0  
#   9    10    VLA:W3    25.0 m   -107.37.08.9  +33.54.00.1  
#   10   11    VLA:N4    25.0 m   -107.37.06.5  +33.54.06.1  
#   11   12    VLA:W5    25.0 m   -107.37.13.0  +33.53.57.8  
#   12   13    VLA:N3    25.0 m   -107.37.06.3  +33.54.04.8  
#   13   14    VLA:N1    25.0 m   -107.37.06.0  +33.54.01.8  
#   14   15    VLA:N2    25.0 m   -107.37.06.2  +33.54.03.5  
#   15   16    VLA:E7    25.0 m   -107.36.52.4  +33.53.56.5  
#   16   17    VLA:E8    25.0 m   -107.36.48.9  +33.53.55.1  
#   17   18    VLA:W4    25.0 m   -107.37.10.8  +33.53.59.1  
#   18   19    VLA:E5    25.0 m   -107.36.58.4  +33.53.58.8  
#   19   20    VLA:W9    25.0 m   -107.37.25.1  +33.53.51.0  
#   20   21    VLA:W6    25.0 m   -107.37.15.6  +33.53.56.4  
#   21   22    VLA:E4    25.0 m   -107.37.00.8  +33.53.59.7  
#   23   24    VLA:E2    25.0 m   -107.37.04.4  +33.54.01.1  
#   24   25    VLA:N6    25.0 m   -107.37.06.9  +33.54.10.3  
#   25   26    VLA:N9    25.0 m   -107.37.07.8  +33.54.19.0  
#   26   27    VLA:N8    25.0 m   -107.37.07.5  +33.54.15.8  
#   27   28    VLA:W7    25.0 m   -107.37.18.4  +33.53.54.8  
#  
# Tables:  
#    MAIN                   22653 rows  
#    ANTENNA                   28 rows  
#    DATA_DESCRIPTION           1 row  
#    DOPPLER             <absent>  
#    FEED                      28 rows  
#    FIELD                      3 rows  
#    FLAG_CMD             <empty>  
#    FREQ_OFFSET         <absent>  
#    HISTORY                  273 rows  
#    OBSERVATION                1 row  
#    POINTING                 168 rows  
#    POLARIZATION               1 row  
#    PROCESSOR            <empty>  
#    SOURCE                     3 rows  
#    SPECTRAL_WINDOW            1 row  
#    STATE                <empty>  
#    SYSCAL              <absent>  
#    WEATHER             <absent>  
#  
#  
#=====================================================================  
#  
# Get rid of the autocorrelations from the MS  
#  
print ’--Flagautocorr--’  
 
# Don’t default this one either, there is only one parameter (vis)  
 
flagautocorr()  
 
#  
#=====================================================================  
#  
# Set the fluxes of the primary calibrator(s)  
#  
print ’--Setjy--’  
default(’setjy’)  
 
vis = msfile  
 
#  
# 1331+305 = 3C286 is our primary calibrator  
# Use the wildcard on the end of the source name  
# since the field names in the MS have inherited the  
# AIPS qualifiers  
field = ’1331+305*’  
 
# This is 1.4GHz D-config and 1331+305 is sufficiently unresolved  
# that we dont need a model image.  For higher frequencies  
# (particularly in A and B config) you would want to use one.  
modimage = ’’  
 
# Setjy knows about this source so we dont need anything more  
 
saveinputs(’setjy’,prefix+’.setjy.saved’)  
 
# Pause script if you are running in scriptmode  
if scriptmode:  
    inp()  
    user_check=raw_input(’Return to continue script\n’)  
 
setjy()  
 
#  
# You should see something like this in the logger and casapy.log file:  
#  
# 1331+30500002_0  spwid=  0  [I=14.76, Q=0, U=0, V=0] Jy, (Perley-Taylor 99)  
#  
# So its using 14.76Jy as the flux of 1331+305 in the single Spectral Window  
# in this MS.  
#  
#=====================================================================  
#  
# Bandpass calibration  
#  
print ’--Bandpass--’  
default(’bandpass’)  
 
# We can first do the bandpass on the single 5min scan on 1331+305  
# At 1.4GHz phase stablility should be sufficient to do this without  
# a first (rough) gain calibration.  This will give us the relative  
# antenna gain as a function of frequency.  
 
vis = msfile  
 
# set the name for the output bandpass caltable  
btable = prefix + ’.bcal’  
caltable = btable  
 
# No gain tables yet  
gaintable = ’’  
gainfield = ’’  
interp = ’’  
 
# Use flux calibrator 1331+305 = 3C286 (FIELD_ID 0) as bandpass calibrator  
field = ’0’  
# all channels  
spw = ’’  
# No other selection  
selectdata = False  
 
# In this band we do not need a-priori corrections for  
# antenna gain-elevation curve or atmospheric opacity  
# (at 8GHz and above you would want these)  
gaincurve = False  
opacity = 0.0  
 
# Choose bandpass solution type  
# Pick standard time-binned B (rather than BPOLY)  
bandtype = ’B’  
 
# set solution interval arbitrarily long (get single bpass)  
solint = ’inf’  
combine = ’scan’  
 
# reference antenna Name 15 (15=VLA:N2) (Id 14)  
refant = ’15’  
 
saveinputs(’bandpass’,prefix+’.bandpass.saved’)  
 
# Pause script if you are running in scriptmode  
if scriptmode:  
    inp()  
    user_check=raw_input(’Return to continue script\n’)  
 
bandpass()  
 
#  
#=====================================================================  
#  
# Use plotcal to examine the bandpass solutions  
#  
print ’--Plotcal (bandpass)--’  
default(’plotcal’)  
 
caltable = btable  
field = ’0’  
 
# Set up 2x1 panels - upper panel amp vs. channel  
subplot = 211  
yaxis = ’amp’  
# No output file yet (wait to plot next panel)  
 
saveinputs(’plotcal’,prefix+’.plotcal.b.amp.saved’)  
 
if scriptmode:  
    showgui = True  
else:  
    showgui = False  
 
plotcal()  
#  
# Set up 2x1 panels - lower panel phase vs. channel  
subplot = 212  
yaxis = ’phase’  
 
saveinputs(’plotcal’,prefix+’.plotcal.b.phase.saved’)  
 
#  
# Note the rolloff in the start and end channels.  Looks like  
# channels 6-56 (out of 0-62) are the best  
 
# Pause script if you are running in scriptmode  
if scriptmode:  
    # If you want to do this interactively and iterate over antenna, set  
    # iteration = ’antenna’  
    showgui = True  
    plotcal()  
    user_check=raw_input(’Return to continue script\n’)  
else:  
    # No GUI for this script  
    showgui = False  
    # Now send final plot to file in PNG format (via .png suffix)  
    figfile = caltable + ’.plotcal.png’  
    plotcal()  
 
#=====================================================================  
#  
# Gain calibration  
#  
print ’--Gaincal--’  
default(’gaincal’)  
 
# Armed with the bandpass, we now solve for the  
# time-dependent antenna gains  
 
vis = msfile  
 
# set the name for the output gain caltable  
gtable = prefix + ’.gcal’  
caltable = gtable  
 
# Use our previously determined bandpass  
# Note this will automatically be applied to all sources  
# not just the one used to determine the bandpass  
gaintable = btable  
gainfield = ’’  
 
# Use nearest (there is only one bandpass entry)  
interp = ’nearest’  
 
# Gain calibrators are 1331+305 and 1445+099 (FIELD_ID 0 and 1)  
field = ’0,1’  
 
# We have only a single spectral window (SPW 0)  
# Choose 51 channels 6-56 out of the 63  
# to avoid end effects.  
# Channel selection is done inside spw  
spw = ’0:6~56’  
 
# No other selection  
selectdata = False  
 
# In this band we do not need a-priori corrections for  
# antenna gain-elevation curve or atmospheric opacity  
# (at 8GHz and above you would want these)  
gaincurve = False  
opacity = 0.0  
 
# scan-based G solutions for both amplitude and phase  
gaintype = ’G’  
solint = ’inf’  
combine = ’’  
calmode = ’ap’  
 
# minimum SNR allowed  
minsnr = 1.0  
 
# reference antenna 15 (15=VLA:N2)  
refant = ’15’  
 
saveinputs(’gaincal’,prefix+’.gaincal.saved’)  
 
if scriptmode:  
    inp()  
    user_check=raw_input(’Return to continue script\n’)  
 
gaincal()  
 
#  
#=====================================================================  
#  
# Bootstrap flux scale  
#  
print ’--Fluxscale--’  
default(’fluxscale’)  
 
vis = msfile  
 
# set the name for the output rescaled caltable  
ftable = prefix + ’.fluxscale’  
fluxtable = ftable  
 
# point to our first gain cal table  
caltable = gtable  
 
# we will be using 1331+305 (the source we did setjy on) as  
# our flux standard reference - note its extended name as in  
# the FIELD table summary above (it has a VLA seq number appended)  
reference = ’1331*’  
 
# we want to transfer the flux to our other gain cal source 1445+099  
transfer = ’1445*’  
 
saveinputs(’fluxscale’,prefix+’.fluxscale.saved’)  
 
# Pause script if you are running in scriptmode  
if scriptmode:  
    inp()  
    user_check=raw_input(’Return to continue script\n’)  
 
fluxscale()  
 
# In the logger you should see something like:  
# Flux density for 1445+09900002_0 in SpW=0 is:  
#     2.48576 +/- 0.00123122 (SNR = 2018.94, nAnt= 27)  
 
# If you run plotcal() on the tablein = ’ngc5921.demo.fluxscale’  
# you will see now it has brought the amplitudes in line between  
# the first scan on 1331+305 and the others on 1445+099  
 
#  
#=====================================================================  
#  
# Now use plotcal to examine the gain solutions  
#  
print ’--Plotcal (fluxscaled gains)--’  
default(’plotcal’)  
 
caltable = ftable  
field = ’0,1’  
 
# Set up 2x1 panels - upper panel amp vs. time  
subplot = 211  
yaxis = ’amp’  
# No output file yet (wait to plot next panel)  
 
saveinputs(’plotcal’,prefix+’.plotcal.gscaled.amp.saved’)  
 
if scriptmode:  
    showgui = True  
else:  
    showgui = False  
 
plotcal()  
#  
# Set up 2x1 panels - lower panel phase vs. time  
subplot = 212  
yaxis = ’phase’  
 
saveinputs(’plotcal’,prefix+’.plotcal.gscaled.phase.saved’)  
 
#  
# The amp and phase coherence looks good  
 
# Pause script if you are running in scriptmode  
if scriptmode:  
    # If you want to do this interactively and iterate over antenna, set  
    #iteration = ’antenna’  
    showgui = True  
    plotcal()  
    user_check=raw_input(’Return to continue script\n’)  
else:  
    # No GUI for this script  
    showgui = False  
    # Now send final plot to file in PNG format (via .png suffix)  
    figfile = caltable + ’.plotcal.png’  
    plotcal()  
 
#=====================================================================  
#  
# Apply our calibration solutions to the data  
# (This will put calibrated data into the CORRECTED_DATA column)  
#  
print ’--ApplyCal--’  
default(’applycal’)  
 
vis = msfile  
 
# We want to correct the calibrators using themselves  
# and transfer from 1445+099 to itself and the target N5921  
 
# Start with the fluxscale/gain and bandpass tables  
gaintable = [ftable,btable]  
 
# pick the 1445+099 out of the gain table for transfer  
# use all of the bandpass table  
gainfield = [’1’,’*’]  
 
# interpolation using linear for gain, nearest for bandpass  
interp = [’linear’,’nearest’]  
 
# only one spw, do not need mapping  
spwmap = []  
 
# all channels  
spw = ’’  
selectdata = False  
 
# as before  
gaincurve = False  
opacity = 0.0  
 
# select the fields for 1445+099 and N5921  
field = ’1,2’  
 
applycal()  
 
# Now for completeness apply 1331+305 to itself  
 
field = ’0’  
gainfield = [’0’,’*’]  
 
# The CORRECTED_DATA column now contains the calibrated visibilities  
 
saveinputs(’applycal’,prefix+’.applycal.saved’)  
 
# Pause script if you are running in scriptmode  
if scriptmode:  
    inp()  
    user_check=raw_input(’Return to continue script\n’)  
 
applycal()  
 
#  
#=====================================================================  
#  
# Now use plotms to plot the calibrated target data (before contsub)  
#  
print ’--Plotms (NGC5921)--’  
default(’plotms’)  
 
vis = msfile  
 
field = ’2’  
# Edge channels are bad  
spw = ’0:4~59’  
 
# Time average across scans  
timebin = ’86000.’  
averagedata=T  
averagescans=T  
 
xaxis = ’channel’  
yaxis = ’amp’  
ydatacolumn = ’corrected’  
# No output file yet (wait to plot next panel)  
 
saveinputs(’plotms’,prefix+’.plotms.final.amp.saved’)  
 
print "Amp averaged across time and baseline (upper)"  
 
plotms()  
#  
yaxis = ’phase’  
ydatacolumn = ’corrected’  
# Time average across scans and baselines  
timebin = ’86000.’  
averagedata=T  
averagescans=T  
averagebaselines=T  
 
saveinputs(’plotms’,prefix+’.plotms.final.phase.saved’)  
 
print "Phase averaged across time and baseline (lower)"  
 
print "Final calibrated data"  
 
 
#=====================================================================  
#  
# Split the sources out, pick off the CORRECTED_DATA column  
#  
#  
# Split NGC5921 data (before continuum subtraction)  
#  
print ’--Split NGC5921 Data--’  
default(’split’)  
 
vis = msfile  
splitms = prefix + ’.src.split.ms’  
outputvis = splitms  
field = ’N5921*’  
spw = ’’  
datacolumn = ’corrected’  
 
saveinputs(’split’,prefix+’.split.n5921.saved’)  
 
split()  
 
print "Created "+splitms  
 
# If you want, split out the calibrater 1445+099 field, all chans  
#print ’--Split 1445+099 Data--’  
#  
#calsplitms = prefix + ’.cal.split.ms’  
#outputvis = calsplitms  
#field = ’1445*’  
#  
#saveinputs(’split’,prefix+’.split.1445.saved’)  
#  
#split()  
 
#=====================================================================  
#  
# Here is how to export the NGC5921 data as UVFITS  
# Start with the split file.  
# Since this is a split dataset, the calibrated data is  
# in the DATA column already.  
# Write as a multisource UVFITS (with SU table)  
# even though it will have only one field in it  
# Run asynchronously so as not to interfere with other tasks  
# (BETA: also avoids crash on next importuvfits)  
#  
#print ’--Export UVFITS--’  
#default(’exportuvfits’)  
#  
#srcuvfits = prefix + ’.split.uvfits’  
#  
#vis = splitms  
#fitsfile = srcuvfits  
#datacolumn = ’data’  
#multisource = True  
#async = True  
#  
#saveinputs(’exportuvfits’,prefix+’.exportuvfits.saved’)  
#  
#myhandle = exportuvfits()  
#  
#print "The return value for this exportuvfits async task for tm is "+str(myhandle)  
 
#=====================================================================  
#  
# UV-plane continuum subtraction on the target  
# use the split ms  
# (this will update the CORRECTED_DATA column)  
#  
print ’--UV Continuum Subtract--’  
default(’uvcontsub’)  
 
vis = splitms  
 
field = ’N5921*’  
# Use channels 4-6 and 50-59 for continuum  
fitspw=’0:4~6;50~59’  
 
# Output all of spw 0  
spw = ’0’  
 
# Averaging time (none)  
solint = 0.0  
 
# Fit only a mean level  
fitorder = 0  
 
# Do the uv-plane subtraction  
fitmode = ’subtract’  
 
# Let it split out the data automatically for us  
splitdata = True  
 
saveinputs(’uvcontsub’,prefix+’.uvcontsub.saved’)  
 
# Pause script if you are running in scriptmode  
if scriptmode:  
    inp()  
    user_check=raw_input(’Return to continue script\n’)  
 
uvcontsub()  
 
# You will see it made two new MS:  
# <vis>.cont  
# <vis>.contsub  
 
srcsplitms = splitms + ’.contsub’  
 
# Note that ngc5921.demo.ms.contsub contains the uv-subtracted  
# visibilities (in its DATA column), and ngc5921.demo.ms.cont  
# the pseudo-continuum visibilities (as fit).  
 
# The original ngc5921.demo.ms now contains the uv-continuum  
# subtracted vis in its CORRECTED_DATA column and the continuum  
# in its MODEL_DATA column as per the fitmode=’subtract’  
 
# Done with calibration  
#=====================================================================  
#  
# Here is how to make a dirty image cube  
#  
#print ’--Clean (dirty image)--’  
#default(’clean’)  
 
# Pick up our split source continuum-subtracted data  
#vis = srcsplitms  
#dirtyname = prefix + ’.dirtyimg’  
#imagename = dirtyname  
#  
#mode = ’channel’  
#nchan = 46  
#start = 5  
#width = 1  
#  
#field = ’0’  
#spw = ’’  
#imsize = [256,256]  
#cell = [15.,15.]  
#weighting = ’briggs’  
#robust = 0.5  
 
# No cleaning  
#niter = 0  
 
#saveinputs(’clean’,prefix+’.invert.saved’)  
 
# Pause script if you are running in scriptmode  
#if scriptmode:  
#    inp()  
#    user_check=raw_input(’Return to continue script\n’)  
#  
#clean()  
 
#dirtyimage = dirtyname+’.image’  
 
# Get the dirty image cube statistics  
#dirtystats = imstat(dirtyimage)  
 
#=====================================================================  
#  
# Now clean an image cube of N5921  
#  
print ’--Clean (clean)--’  
default(’clean’)  
 
# Pick up our split source continuum-subtracted data  
vis = srcsplitms  
 
# Make an image root file name  
imname = prefix + ’.cleanimg’  
imagename = imname  
 
# Set up the output image cube  
mode = ’channel’  
nchan = 46  
start = 5  
width = 1  
 
# This is a single-source MS with one spw  
field = ’0’  
spw = ’’  
 
# Standard gain factor 0.1  
gain = 0.1  
 
# Set the output image size and cell size (arcsec)  
imsize = [256,256]  
 
# Do a simple Clark clean  
psfmode = ’clark’  
# No Cotton-Schwab iterations  
csclean = False  
 
# If desired, you can do a Cotton-Schwab clean  
# but will have only marginal improvement for this data  
#csclean = True  
# Twice as big for Cotton-Schwab (cleans inner quarter)  
#imsize = [512,512]  
 
# Pixel size 15 arcsec for this data (1/3 of 45" beam)  
# VLA D-config L-band  
cell = [15.,15.]  
 
# Fix maximum number of iterations  
niter = 6000  
 
# Also set flux residual threshold (in mJy)  
threshold=8.0  
 
# Set up the weighting  
# Use Briggs weighting (a moderate value, on the uniform side)  
weighting = ’briggs’  
robust = 0.5  
 
# Set a cleanbox +/-20 pixels around the center 128,128  
mask = [108,108,148,148]  
 
# But if you had a cleanbox saved in a file, e.g. "regionfile.txt"  
# you could use it:  
#mask=’regionfile.txt’  
#  
# If you don’t want any clean boxes or masks, then  
#mask = ’’  
 
# If you want interactive clean set to True  
#interactive=True  
interactive=False  
 
saveinputs(’clean’,prefix+’.clean.saved’)  
 
# Pause script if you are running in scriptmode  
if scriptmode:  
    inp()  
    user_check=raw_input(’Return to continue script\n’)  
 
clean()  
 
# Should find stuff in the logger like:  
#  
# Fitted beam used in restoration: 51.5643 by 45.6021 (arcsec)  
#     at pa 14.5411 (deg)  
#  
# It will have made the images:  
# -----------------------------  
# ngc5921.demo.cleanimg.flux  
# ngc5921.demo.cleanimg.image  
# ngc5921.demo.cleanimg.mask  
# ngc5921.demo.cleanimg.model  
# ngc5921.demo.cleanimg.psf  
# ngc5921.demo.cleanimg.residual  
 
clnimage = imname+’.image’  
 
#=====================================================================  
#  
# Done with imaging  
# Now view the image cube of N5921  
#  
if scriptmode:  
    print ’--View image--’  
    print "Use Spectral Profile Tool to get line profile in box in center"  
    viewer(clnimage)  
    user_check=raw_input(’Return to continue script\n’)  
 
#=====================================================================  
#  
# Here is how to export the Final CLEAN Image as FITS  
# Run asynchronously so as not to interfere with other tasks  
# (BETA: also avoids crash on next importfits)  
#  
#print ’--Final Export CLEAN FITS--’  
#default(’exportfits’)  
#  
#clnfits = prefix + ’.cleanimg.fits’  
#  
#imagename = clnimage  
#fitsimage = clnfits  
#async = True  
#  
#saveinputs(’exportfits’,prefix+’.exportfits.saved’)  
#  
#myhandle2 = exportfits()  
#  
#print "The return value for this exportfits async task for tm is "+str(myhandle2)  
 
#=====================================================================  
#  
# Print the image header  
#  
print ’--Imhead--’  
default(’imhead’)  
 
imagename = clnimage  
 
mode = ’summary’  
 
imhead()  
 
# A summary of the cube will be seen in the logger  
 
#=====================================================================  
#  
# Get the cube statistics  
#  
print ’--Imstat (cube)--’  
default(’imstat’)  
 
imagename = clnimage  
 
# Do whole image  
box = ’’  
# or you could stick to the cleanbox  
#box = ’108,108,148,148’  
 
cubestats = imstat()  
 
# Statistics will printed to the terminal, and the output  
# parameter will contain a dictionary of the statistics  
 
#=====================================================================  
#  
# Get some image moments  
#  
print ’--ImMoments--’  
default(’immoments’)  
 
imagename = clnimage  
 
# Do first and second moments  
moments = [0,1]  
 
# Need to mask out noisy pixels, currently done  
# using hard global limits  
excludepix = [-100,0.009]  
 
# Collapse along the spectral (channel) axis  
axis  = ’spectral’  
# Include all planes  
chans = ’’  
 
# Output root name  
momfile = prefix + ’.moments’  
outfile = momfile  
 
saveinputs(’immoments’,prefix+’.immoments.saved’)  
 
# Pause script if you are running in scriptmode  
if scriptmode:  
    inp()  
    user_check=raw_input(’Return to continue script\n’)  
 
immoments()  
 
momzeroimage = momfile + ’.integrated’  
momoneimage = momfile + ’.weighted_coord’  
 
# It will have made the images:  
# --------------------------------------  
# ngc5921.demo.moments.integrated  
# ngc5921.demo.moments.weighted_coord  
 
#  
#=====================================================================  
#  
# Get some statistics of the moment images  
#  
print ’--Imstat (moments)--’  
default(’imstat’)  
 
imagename = momzeroimage  
momzerostats = imstat()  
 
imagename = momoneimage  
momonestats = imstat()  
 
#=====================================================================  
#  
# Now view the moments  
#  
if scriptmode:  
    print ’--View image (Moments)--’  
    viewer(momzeroimage)  
    print "You can add mom-1 image "+momoneimage  
    print "as a contour plot"  
    user_check=raw_input(’Return to continue script\n’)  
 
#=====================================================================  
#  
# Set up an output logfile  
import datetime  
datestring=datetime.datetime.isoformat(datetime.datetime.today())  
 
outfile = ’out.’+prefix+’.’+datestring+’.log’  
logfile=open(outfile,’w’)  
print >>logfile,’Results for ’+prefix+’ :’  
print >>logfile,""  
 
#=====================================================================  
#  
# Can do some image statistics if you wish  
# Treat this like a regression script  
# WARNING: currently requires toolkit  
#  
print ’ NGC5921 results ’  
print ’ =============== ’  
 
print >>logfile,’ NGC5921 results ’  
print >>logfile,’ =============== ’  
 
#  
# Use the ms tool to get max of the MSs  
# Eventually should be available from a task  
#  
# Pull the max cal amp value out of the MS (if you split this)  
#ms.open(calsplitms)  
#thistest_cal = max(ms.range(["amplitude"]).get(’amplitude’))  
#ms.close()  
#oldtest_cal = 34.0338668823  
#diff_cal = abs((oldtest_cal-thistest_cal)/oldtest_cal)  
#  
#print ’ Calibrator data ampl max = ’,thistest_cal  
#print ’   Previous: cal data max = ’,oldtest_cal  
#print ’   Difference (fractional) = ’,diff_cal  
#print ’’  
#  
#print >>logfile,’ Calibrator data ampl max = ’,thistest_cal  
#print >>logfile,’   Previous: cal data max = ’,oldtest_cal  
#print >>logfile,’   Difference (fractional) = ’,diff_cal  
#print >>logfile,’’  
 
# Pull the max src amp value out of the MS  
ms.open(srcsplitms)  
thistest_src = max(ms.range(["amplitude"]).get(’amplitude’))  
ms.close()  
oldtest_src =  46.2060050964 # now in all chans  
diff_src = abs((oldtest_src-thistest_src)/oldtest_src)  
 
print ’ Target Src data ampl max = ’,thistest_src  
print ’   Previous: src data max = ’,oldtest_src  
print ’   Difference (fractional) = ’,diff_src  
print ’’  
 
print >>logfile,’ Target Src data ampl max = ’,thistest_src  
print >>logfile,’   Previous: src data max = ’,oldtest_src  
print >>logfile,’   Difference (fractional) = ’,diff_src  
print >>logfile,’’  
 
#  
# Now use the stats produced by imstat above  
#  
# DIRTY IMAGE MAX & RMS (IF YOU MADE A DIRTY IMAGE)  
#  
#thistest_dirtymax=dirtystats[’max’][0]  
#oldtest_dirtymax = 0.0515365377069  
#diff_dirtymax = abs((oldtest_dirtymax-thistest_dirtymax)/oldtest_dirtymax)  
#  
#print ’ Dirty image max = ’,thistest_dirtymax  
#print ’   Previous: max = ’,oldtest_dirtymax  
#print ’   Difference (fractional) = ’,diff_dirtymax  
#print ’’  
#  
#print >>logfile,’ Dirty Image max = ’,thistest_dirtymax  
#print >>logfile,’   Previous: max = ’,oldtest_dirtymax  
#print >>logfile,’   Difference (fractional) = ’,diff_dirtymax  
#print >>logfile,’’  
#  
#thistest_dirtyrms=dirtystats[’rms’][0]  
#oldtest_dirtyrms = 0.00243866862729  
#diff_dirtyrms = abs((oldtest_dirtyrms-thistest_dirtyrms)/oldtest_dirtyrms)  
#  
#print ’ Dirty image rms = ’,thistest_dirtyrms  
#print ’   Previous: rms = ’,oldtest_dirtyrms  
#print ’   Difference (fractional) = ’,diff_dirtyrms  
#print ’’  
#  
#print >>logfile,’ Dirty Image rms = ’,thistest_dirtyrms  
#print >>logfile,’   Previous: rms = ’,oldtest_dirtyrms  
#print >>logfile,’   Difference (fractional) = ’,diff_dirtyrms  
#print >>logfile,’’  
 
# Now the clean image  
#  
thistest_immax=cubestats[’max’][0]  
oldtest_immax = 0.052414759993553162  
diff_immax = abs((oldtest_immax-thistest_immax)/oldtest_immax)  
 
print ’ Clean image max = ’,thistest_immax  
print ’   Previous: max = ’,oldtest_immax  
print ’   Difference (fractional) = ’,diff_immax  
print ’’  
 
print >>logfile,’ Clean Image max = ’,thistest_immax  
print >>logfile,’   Previous: max = ’,oldtest_immax  
print >>logfile,’   Difference (fractional) = ’,diff_immax  
print >>logfile,’’  
 
thistest_imrms=cubestats[’rms’][0]  
oldtest_imrms = 0.0020218724384903908  
diff_imrms = abs((oldtest_imrms-thistest_imrms)/oldtest_imrms)  
 
print ’ Clean image rms = ’,thistest_imrms  
print ’   Previous: rms = ’,oldtest_imrms  
print ’   Difference (fractional) = ’,diff_imrms  
print ’’  
 
print >>logfile,’ Clean image rms = ’,thistest_imrms  
print >>logfile,’   Previous: rms = ’,oldtest_imrms  
print >>logfile,’   Difference (fractional) = ’,diff_imrms  
print >>logfile,’’  
 
# Now the moment images  
#  
thistest_momzeromax=momzerostats[’max’][0]  
oldtest_momzeromax = 1.40223777294  
diff_momzeromax = abs((oldtest_momzeromax-thistest_momzeromax)/oldtest_momzeromax)  
 
print ’ Moment 0 image max = ’,thistest_momzeromax  
print ’   Previous: m0 max = ’,oldtest_momzeromax  
print ’   Difference (fractional) = ’,diff_momzeromax  
print ’’  
 
print >>logfile,’ Moment 0 image max = ’,thistest_momzeromax  
print >>logfile,’   Previous: m0 max = ’,oldtest_momzeromax  
print >>logfile,’   Difference (fractional) = ’,diff_momzeromax  
print >>logfile,’’  
 
thistest_momoneavg=momonestats[’mean’][0]  
oldtest_momoneavg = 1479.77119646  
diff_momoneavg = abs((oldtest_momoneavg-thistest_momoneavg)/oldtest_momoneavg)  
 
print ’ Moment 1 image mean = ’,thistest_momoneavg  
print ’   Previous: m1 mean = ’,oldtest_momoneavg  
print ’   Difference (fractional) = ’,diff_momoneavg  
print ’’  
print ’--- Done ---’  
 
print >>logfile,’ Moment 1 image mean = ’,thistest_momoneavg  
print >>logfile,’   Previous: m1 mean = ’,oldtest_momoneavg  
print >>logfile,’   Difference (fractional) = ’,diff_momoneavg  
print >>logfile,’’  
print >>logfile,’--- Done ---’  
 
# Should see output like:  
#  
# Clean image max should be  0.0524147599936  
# Found : Image Max =  0.0523551553488  
# Difference (fractional) =  0.00113717290288  
#  
# Clean image rms should be  0.00202187243849  
# Found : Image rms =  0.00202226242982  
# Difference (fractional) =  0.00019288621809  
#  
# Moment 0 image max should be  1.40223777294  
# Found : Moment 0 Max =  1.40230333805  
# Difference (fractional) =  4.67574844349e-05  
#  
# Moment 1 image mean should be  1479.77119646  
# Found : Moment 1 Mean =  1479.66974528  
# Difference (fractional) =  6.85586935973e-05  
#  
#=====================================================================  
# Done  
#  
logfile.close()  
print "Results are in "+outfile