######################################################################
# #
# 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 #
# Last Updated STM 2007-10-10 (Beta) add export #
# #
######################################################################
import time
import os
#
# Set up some useful variables
#
# 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'
# The prefix to use for all output files
prefix='ngc5921.usecase'
# Clean up old 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
importuvfits()
#
# Note that there will be a ngc5921.usecase.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 casa.log file
# something like:
#
# MeasurementSet Name: /home/sandrock2/smyers/Testing2/Sep07/ngc5921.usecase.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 '--Flag auto-correlation--'
default('flagdata')
vis = msfile
mode = 'manual'
autocorr = True
flagdata()
#
#=====================================================================
#
# 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
setjy()
#
# You should see something like this in the logger and casa.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
# Choose bandpass solution type
# Pick standard time-binned B (rather than BPOLY)
bandtype = 'B'
# set solution interval arbitrarily long (get single bpass)
solint = 86400.0
# reference antenna Name 15 (15=VLA:N2) (Id 14)
refant = '15'
bandpass()
# You can use plotcal to examine the solutions
#default('plotcal')
#tablein = btable
#yaxis = 'amp'
#field = '0'
#multiplot = True
#plotcal()
#
#yaxis = 'phase'
#plotcal()
#
# Note the rolloff in the start and end channels. Looks like
# channels 6-56 (out of 0-62) are the best
#=====================================================================
#
# 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
# scan-based G solutions for both amplitude and phase
gaintype = 'G'
solint = 0.
calmode = 'ap'
# minimum SNR allowed
minsnr = 1.0
# reference antenna 15 (15=VLA:N2)
refant = '15'
gaincal()
# You can use plotcal to examine the gain solutions
#default('plotcal')
#tablein = gtable
#yaxis = 'amp'
#field = '0,1'
#multiplot = True
#plotcal()
#
#yaxis = 'phase'
#plotcal()
#
# The amp and phase coherence looks good
#=====================================================================
#
# 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*'
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.usecase.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
#=====================================================================
#
# 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
# 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
applycal()
#=====================================================================
#
# Split the gain calibrater data, then the target
#
print '--Split 1445+099 Data--'
default('split')
vis = msfile
# We first want to write out the corrected data for the calibrator
# Make an output vis file
calsplitms = prefix + '.cal.split.ms'
outputvis = calsplitms
# Select the 1445+099 field, all chans
field = '1445*'
spw = ''
# pick off the CORRECTED_DATA column
datacolumn = 'corrected'
split()
#
# Now split NGC5921 data (before continuum subtraction)
#
print '--Split NGC5921 Data--'
splitms = prefix + '.src.split.ms'
outputvis = splitms
# Pick off N5921
field = 'N5921*'
split()
#=====================================================================
#
# Export the NGC5921 data as UVFITS
# Start with the split file.
#
print '--Export UVFITS--'
default('exportuvfits')
srcuvfits = prefix + '.split.uvfits'
vis = splitms
fitsfile = srcuvfits
# Since this is a split dataset, the calibrated data is
# in the DATA column already.
datacolumn = 'data'
# Write as a multisource UVFITS (with SU table)
# even though it will have only one field in it
multisource = True
exportuvfits()
#=====================================================================
#
# UV-plane continuum subtraction on the target
# (this will update the CORRECTED_DATA column)
#
print '--UV Continuum Subtract--'
default('uvcontsub')
vis = msfile
# Pick off N5921
field = 'N5921*'
# Use channels 4-6 and 50-59 for continuum
#spw = '0:4~6;50~59'
# BETA ALERT: still does not use standard notation
spw = '0'
channels = range(4,7)+range(50,60)
# 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
uvcontsub()
# You will see it made two new MS:
# ngc5921.usecase.ms.cont
# ngc5921.usecase.ms.contsub
srcsplitms = msfile + '.contsub'
# Note that ngc5921.usecase.ms.contsub contains the uv-subtracted
# visibilities (in its DATA column), and ngc5921.usecase.ms.contsub
# the pseudo-continuum visibilities (as fit).
# The original ngc5921.usecase.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
# Now clean an image cube of N5921
#
print '--Clean--'
default('clean')
# Pick up our split source data
vis = srcsplitms
# Make an image root file name
imname = prefix + '.clean'
imagename = imname
# Set up the output image cube
mode = 'channel'
nchan = 46
start = 5
step = 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
alg = 'clark'
# If desired, you can do a Cotton-Schwab clean
# but will have only marginal improvement for this data
#alg = 'csclean'
# 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'
rmode = 'norm'
robust = 0.5
# No clean mask or cleanbox for now
mask = ''
cleanbox = []
# But if you had a cleanbox saved in a file, e.g. "regionfile.txt"
# you could use it:
#cleanbox='regionfile.txt'
#
# and if you wanted to use interactive clean
#cleanbox='interactive'
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.usecase.clean.image
# ngc5921.usecase.clean.model
# ngc5921.usecase.clean.residual
clnimage = imname+'.image'
#=====================================================================
#
# Done with imaging
# Now view the image cube of N5921
#
#print '--View image--'
#viewer(clnimage,'image')
#=====================================================================
#
# Export the Final CLEAN Image as FITS
#
print '--Final Export CLEAN FITS--'
default('exportfits')
clnfits = prefix + '.clean.fits'
imagename = clnimage
fitsimage = clnfits
exportfits()
#=====================================================================
#
# Get some image statistics
#
print '--Imhead--'
default('imhead')
imagename = clnimage
mode = 'stats'
cubestats = imhead()
# A summary of the cube will be seen in the logger
# cubestats 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]
# Include all planes
planes = ''
# Output root name
momfile = prefix + '.moments'
outfile = momfile
immoments()
momzeroimage = momfile + '.integrated'
momoneimage = momfile + '.weighted_coord'
# It will have made the images:
# --------------------------------------
# ngc5921.usecase.moments.integrated
# ngc5921.usecase.moments.weighted_coord
#=====================================================================
#
# Get some statistics of the moment images
#
print '--Imhead--'
default('imhead')
mode = 'stats'
imagename = momzeroimage
momzerostats = imhead()
imagename = momoneimage
momonestats = imhead()
#=====================================================================
#
# Can do some image statistics if you wish
# Treat this like a regression script
# WARNING: currently requires toolkit
#
print ' NGC5921 results '
print ' =============== '
#
# 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
ms.open(calsplitms)
thistest_cal = max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
oldtest_cal = 34.0338668823
print ' Cal Max amplitude should be ',oldtest_cal
print ' Found : Max = ',thistest_cal
diff_cal = abs((oldtest_cal-thistest_cal)/oldtest_cal)
print ' Difference (fractional) = ',diff_cal
print ''
# 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 = 1.37963354588 # This was in chans 5-50
oldtest_src = 46.2060050964 # now in all chans
print ' Src Max amplitude should be ',oldtest_src
print ' Found : Max = ',thistest_src
diff_src = abs((oldtest_src-thistest_src)/oldtest_src)
print ' Difference (fractional) = ',diff_src
#
# Now use the stats produced by imhead above
#
print ''
# Pull the max from the cubestats dictionary
# created above using imhead
thistest_immax=cubestats['max'][0]
oldtest_immax = 0.052414759993553162
print ' Clean image max should be ',oldtest_immax
print ' Found : Image Max = ',thistest_immax
diff_immax = abs((oldtest_immax-thistest_immax)/oldtest_immax)
print ' Difference (fractional) = ',diff_immax
print ''
# Pull the rms from the cubestats dictionary
thistest_imrms=cubestats['rms'][0]
oldtest_imrms = 0.0020218724384903908
print ' Clean image rms should be ',oldtest_imrms
print ' Found : Image rms = ',thistest_imrms
diff_imrms = abs((oldtest_imrms-thistest_imrms)/oldtest_imrms)
print ' Difference (fractional) = ',diff_imrms
print ''
# Pull the max from the momzerostats dictionary
thistest_momzeromax=momzerostats['max'][0]
oldtest_momzeromax = 1.40223777294
print ' Moment 0 image max should be ',oldtest_momzeromax
print ' Found : Moment 0 Max = ',thistest_momzeromax
diff_momzeromax = abs((oldtest_momzeromax-thistest_momzeromax)/oldtest_momzeromax)
print ' Difference (fractional) = ',diff_momzeromax
print ''
# Pull the mean from the momonestats dictionary
thistest_momoneavg=momonestats['mean'][0]
oldtest_momoneavg = 1479.77119646
print ' Moment 1 image mean should be ',oldtest_momoneavg
print ' Found : Moment 1 Mean = ',thistest_momoneavg
diff_momoneavg = abs((oldtest_momoneavg-thistest_momoneavg)/oldtest_momoneavg)
print ' Difference (fractional) = ',diff_momoneavg
print ''
print '--- Done ---'
#
#=====================================================================