#-- Self-calibration on continuum, applied to all channels/spw.
#-- Imaging of continuum, continuum subtraction, make cubes for
#-- water maser and another line in the other spw.

#-- Test channel 1:636 containing extended emission is imaged after
#-- continuum subtraction, to compare with other self-cal routes

#-- You could try using the contours of the final continuum image
#-- (made with all 3 EBs) as a guide for the first mask.

#-- Set this option to False if you do not want to do interactive clean.
#-- Note that better results will be obtained if you do make masks interatively,
#-- expanding the area as fainter emission emerges
cleaninteractive=False

#-- Set this option to true if you want to make diagnostic plots along the way.
#-- Note that this will slow down the reduction.
#----- Default is no plots (you can make them later anyway).
calplots=True

refantenna='DV15'

# Self-cal on continuum

# This selection was made after calibration on the maser and imaging but 
# you could make a selection from the visibility spectrum before self-cal
contchans='0:5~224;365~379;425~479;615~739;1185~1249;1460~1569;1615~1699;1820~1919,1:0~214;440~494;580~609;750~789;830~1049;1305~1364;1490~1529;1610~1659;1775~1829;1815~1834;1895~1919'

thesteps=[]
step_title = {0: 'List the data set',
              1: 'Plot the visibility spectrum and peak; select continuum',
              2: 'Make first image of continuum',
              3: 'Self-calibrate cycle 1, phase-only, re-image',
              4: 'Self-calibrate cycle 2, phase-only, re-image',
              5: 'Self-calibrate cycle 3, amp & phase, re-image',
              6: 'Self-calibrate cycle 4, amp & phase, re-image',
              7: 'Make low-res cubes to identify continuum \n(this step can be omitted if you trust our channel selection)',
              8: 'Image continuum',
              9: 'Subtract continuum and make image of spw 1 test channel',
              10: 'Image a line in spw 0',
              11: 'Image a line in spw 1',
}
 
# The Python variable 'mysteps' will control which steps
# are executed when you start the script using
#   execfile('scriptForCalibration.py')
# e.g. setting
#   mysteps = [2,3,4]# before starting the script will make the script execute
# only steps 2, 3, and 4
# Setting mysteps = [] will make it execute all steps.

try:
  print 'List of steps to be executed ...', mysteps
  thesteps = mysteps
except:
  print 'global variable mysteps not set.'
if (thesteps==[]):
  thesteps = range(0,len(step_title))
  print 'Executing all steps: ', thesteps

"""
# PREPARATION
# Select just one of the 3 EBs, with Tsys, WVR, bandpass and phase-ref
# solutions applied.  Shift to constant LSR velocity and average target

basename=['uid___A002_X6d0f96_X1de2']
refantenna='DV15'

vislist=[]
for name in basename:
    for s in range(2):  
        os.system('rm -rf '+name+'spw'+str(s)+'_VYCMa_325.ms')
        mstransform(vis=name+'.ms.split',
        spw=str(s),
        outputvis=name+'spw'+str(s)+'_VYCMa_325.ms',datacolumn='corrected',
        field='V*',
        regridms=True,
        width=2,
        outframe='lsrk',
        timeaverage=True,
        timebin='12.1s')
        vislist.append(name+'spw'+str(s)+'_VYCMa_325.ms')

concat(vis=vislist,
       concatvis='X1de2_VYCMa_325.ms')

"""
mystep = 0 
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# Copy ms and clear old calibration if necessary
  os.system('cp -r X1de2_VYCMa_325.ms X1de2_VYCMa_contcal_325.ms')
  delmod(vis='X1de2_VYCMa_contcal_325.ms', scr=True)
  clearcal(vis='X1de2_VYCMa_contcal_325.ms')

  os.system('rm X1de2_VYCMa_contcal_325_listobs.txt')
  listobs(vis='X1de2_VYCMa_contcal_325.ms', listfile='X1de2_VYCMa_contcal_325_listobs.txt', verbose=True)
  print '<< Printed listobs output to X1de2_VYCMa_contcal_325_listobs.txt'

mystep = 1 
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  for s in ['0','1']:
      os.system('rm -rf X1de2_VYCMa_contcal_325_spw'+s+'_1st-vis-spectrum.png')
      plotms(vis='X1de2_VYCMa_contcal_325.ms', 
             xaxis='frequency', yaxis='amp',
             selectdata=True, spw=s, 
             avgtime='1e8',avgscan=True,avgbaseline=True,
             coloraxis='baseline',
             showgui=False,
             plotfile='X1de2_VYCMa_contcal_325_spw'+s+'_1st-vis-spectrum.png')

# You will notice the very bright maser line in spw 0:981
# even for continuum it is worth plotting this to reveal bad data

  os.system('rm -rf X1de2_VYCMa_contcal_325_0-981_1st_amp.png')
  plotms(vis='X1de2_VYCMa_contcal_325.ms', 
         xaxis='time', yaxis='amp',
         selectdata=True, 
         antenna='DV15&*',
         spw='0:981',
         coloraxis='baseline',
         showgui=False,
         plotfile='X1de2_VYCMa_contcal_325_0-981_1st_amp.png')

  os.system('rm -rf X1de2_VYCMa_contcal_325_0-981_1st_phase.png')
  plotms(vis='X1de2_VYCMa_contcal_325.ms', 
         xaxis='time', yaxis='phase',
         selectdata=True, 
         antenna='DV15&*',
         spw='0:981',
         coloraxis='baseline',
         showgui=False,
         plotfile='X1de2_VYCMa_contcal_325_0-981_1st_phase.png')

# * DA50  has fast phase but it turns out it can be corrected
######################################################
#  Make first image of continuum, omitting dubious data
mystep = 2
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf X1de2_VYCMa_contcal_325_cont0.clean*')
  clean(vis = 'X1de2_VYCMa_contcal_325.ms',
        imagename='X1de2_VYCMa_contcal_325_cont0.clean',
        field='Vy Cma',
        antenna='!DA50', 
        spw=contchans,
        mode='mfs',
        mask=['ellipse[[250pix,280pix],[70pix,67pix],90deg]','ellipse[[270pix,310pix],[95pix,42pix],20deg]'],
        cell='0.01arcsec',npercycle=10,imsize=512,
        psfmode='hogbom',interactive=cleaninteractive)


# Get image statistics
  imin='X1de2_VYCMa_contcal_325_cont0.clean.image'
  noise=imstat(imagename=imin,
             box='250,50,500,200')['rms'][0]

  peak=imstat(imagename=imin,
            box='200,200,330,330')['max'][0]

  print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)   
# rms 0.003 , peak 0.174 , snr 5    

######################################################
#  Self-calibrate, phase-only first, using all data
mystep = 3
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# ft model into data otherwise it seems only one spw gets calibrated
  ft(vis = 'X1de2_VYCMa_contcal_325.ms',
     model='X1de2_VYCMa_contcal_325_cont0.clean.model',
     usescratch=True)

# decrease minsnr to avoid failing solutions - if you shorten solint
# DA50 is not properly corrected.

  os.system('rm -rf X1de2_VYCMa_contcal_325_cont*.ph1')
  gaincal(vis = 'X1de2_VYCMa_contcal_325.ms',
          field='Vy Cma',
          refant=refantenna,
          caltable='X1de2_VYCMa_contcal_325_cont.ph1',
          spw=contchans,
          calmode='p',
          minsnr=2,
          solint='60s')


  if (calplots):
      os.system('rm -rf X1de2_VYCMa_contcal_325_cont*.ph1*.png')
      for ant in range(3):
          plotcal(caltable='X1de2_VYCMa_contcal_325_cont.ph1',
              xaxis = 'time', yaxis = 'phase',
              antenna=str(ant*8)+'~'+str(ant*8+7),    
              subplot=421,
              fontsize=7.0,
              iteration='antenna',
              figfile='X1de2_VYCMa_contcal_325_cont.ph1_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png')

  applycal(vis = 'X1de2_VYCMa_contcal_325.ms',
           field='Vy Cma',
           spw='0,1',
           spwmap=[0,1],
           gaintable='X1de2_VYCMa_contcal_325_cont.ph1',
           calwt = False,
           applymode='calonly',
           flagbackup = False)

  os.system('rm -rf X1de2_VYCMa_contcal_325_cont_ph1.clean*')
  clean(vis = 'X1de2_VYCMa_contcal_325.ms',
        imagename='X1de2_VYCMa_contcal_325_cont_ph1.clean',
        field='Vy Cma', 
        spw=contchans,
        mode='mfs',
        cell='0.01arcsec',npercycle=50,imsize=512,
        psfmode='hogbom',
        mask=['ellipse[[250pix,280pix],[70pix,67pix],90deg]','ellipse[[270pix,310pix],[95pix,42pix],20deg]'],
        interactive=cleaninteractive,
        niter=250)

# Get image statistics
  imin='X1de2_VYCMa_contcal_325_cont_ph1.clean.image'
  noise=imstat(imagename=imin,
               box='250,50,500,200')['rms'][0]

  peak=imstat(imagename=imin,
              box='200,200,330,330')['max'][0]

  print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)   
#rms    0.002, peak  0.197, snr   101


######################################################
#  Self-calibrate cycle 2, phase-only, shorter solint, using all data
mystep = 4
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# ft model into data otherwise it seems only one spw gets calibrated
  ft(vis = 'X1de2_VYCMa_contcal_325.ms',
     model='X1de2_VYCMa_contcal_325_cont_ph1.clean.model',
     usescratch=True)

  os.system('rm -rf X1de2_VYCMa_contcal_325_cont.ph2')
  gaincal(vis = 'X1de2_VYCMa_contcal_325.ms',
          field='Vy Cma',
          refant=refantenna,
          caltable='X1de2_VYCMa_contcal_325_cont.ph2',
          spw=contchans,
          calmode='p',
          solint='36s')

  if (calplots):
      os.system('rm -rf X1de2_VYCMa_contcal_325_cont.ph2*.png')
      for ant in range(3):
          plotcal(caltable='X1de2_VYCMa_contcal_325_cont.ph2',
              xaxis = 'time', yaxis = 'phase',
              antenna=str(ant*8)+'~'+str(ant*8+7),    
              subplot=421,
              fontsize=7.0,
              iteration='antenna',
              figfile='X1de2_VYCMa_contcal_325_cont.ph2_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png')

  applycal(vis = 'X1de2_VYCMa_contcal_325.ms',
           field='Vy Cma',
           spw='0,1',
           spwmap=[0,1],
           gaintable='X1de2_VYCMa_contcal_325_cont.ph2',
           calwt = False,
           applymode='calonly',
           flagbackup = False)

  os.system('rm -rf X1de2_VYCMa_contcal_325_cont_ph2.clean*')
  clean(vis = 'X1de2_VYCMa_contcal_325.ms',
        imagename='X1de2_VYCMa_contcal_325_cont_ph2.clean',
        field='Vy Cma', 
        spw=contchans,
        mode='mfs',
        cell='0.01arcsec',npercycle=50,imsize=512,
        psfmode='hogbom',
        mask=['ellipse[[250pix,280pix],[70pix,67pix],90deg]','ellipse[[270pix,310pix],[95pix,42pix],20deg]'],
        interactive=cleaninteractive,
        niter=400)

# Get image statistics
  imin='X1de2_VYCMa_contcal_325_cont_ph2.clean.image'
  noise=imstat(imagename=imin,
               box='250,50,500,200')['rms'][0]

  peak=imstat(imagename=imin,
              box='200,200,330,330')['max'][0]

  print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)   
# rms    0.002, peak  0.203, snr 111


######################################################
#  Self-calibrate cycle 3, a&p, using all data
mystep = 5
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  ft(vis = 'X1de2_VYCMa_contcal_325.ms',
     model='X1de2_VYCMa_contcal_325_cont_ph2.clean.model',
     usescratch=True)

  os.system('rm -rf X1de2_VYCMa_contcal_325_cont.ap1')
  gaincal(vis = 'X1de2_VYCMa_contcal_325.ms',
          field='Vy Cma',
          refant=refantenna,
          caltable='X1de2_VYCMa_contcal_325_cont.ap1',
          spw=contchans,
          gaintable='X1de2_VYCMa_contcal_325_cont.ph2',
          spwmap=[0,1],
          calmode='ap',
          solint='120s')

  if (calplots):
      os.system('rm -rf X1de2_VYCMa_contcal_325_cont.ap1*.png')
      for ant in range(3):
          plotcal(caltable='X1de2_VYCMa_contcal_325_cont.ap1',
              xaxis = 'time', yaxis = 'amp',
              antenna=str(ant*8)+'~'+str(ant*8+7),    
              subplot=421,
              fontsize=7.0,
              iteration='antenna',
              figfile='X1de2_VYCMa_contcal_325_cont.ap1_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png')

  applycal(vis = 'X1de2_VYCMa_contcal_325.ms',
           field='Vy Cma',
           spw='0,1',
           spwmap=[[0,1],[0,1]],
           applymode='calonly',
           gaintable=['X1de2_VYCMa_contcal_325_cont.ph2','X1de2_VYCMa_contcal_325_cont.ap1'],
           calwt = False,
           flagbackup = False)

#
  os.system('rm -rf X1de2_VYCMa_contcal_325_cont_ap1.clean*')
  clean(vis = 'X1de2_VYCMa_contcal_325.ms',
        imagename='X1de2_VYCMa_contcal_325_cont_ap1.clean',
        field='Vy Cma', 
        spw=contchans,
        mode='mfs',
        cell='0.01arcsec',npercycle=50,imsize=512,
        psfmode='hogbom',
        mask=['ellipse[[250pix,280pix],[70pix,67pix],90deg]','ellipse[[270pix,310pix],[95pix,42pix],20deg]'],
        interactive=cleaninteractive,
        niter=400)

# Get image statistics
  imin='X1de2_VYCMa_contcal_325_cont_ap1.clean.image'
  noise=imstat(imagename=imin,
               box='250,50,500,200')['rms'][0]

  peak=imstat(imagename=imin,
              box='200,200,330,330')['max'][0]

  print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)   
#   rms    0.001, peak    0.217, snr   187


######################################################
#  Self-calibrate cycle 4, ap, short solint, using all data
mystep = 6
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  ft(vis = 'X1de2_VYCMa_contcal_325.ms',
     model='X1de2_VYCMa_contcal_325_cont_ap1.clean.model',
     usescratch=True)

  os.system('rm -rf X1de2_VYCMa_contcal_325_cont.ap2')
  gaincal(vis = 'X1de2_VYCMa_contcal_325.ms',
          field='Vy Cma',
          refant=refantenna,
          caltable='X1de2_VYCMa_contcal_325_cont.ap2',
          spw=contchans,
          gaintable='X1de2_VYCMa_contcal_325_cont.ph2',
          spwmap=[0,1],
          calmode='ap',
          gaintype='T',            #don't fail solutions
          solint='60s')

  if (calplots):
      os.system('rm -rf X1de2_VYCMa_contcal_325_cont.ap2*.png')
      for ant in range(3):
          plotcal(caltable='X1de2_VYCMa_contcal_325_cont.ap2',
              xaxis = 'time', yaxis = 'amp',
              antenna=str(ant*8)+'~'+str(ant*8+7),    
              subplot=421,
              fontsize=7.0,
              iteration='antenna',
              figfile='X1de2_VYCMa_contcal_325_cont.ap2_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png')

  applycal(vis = 'X1de2_VYCMa_contcal_325.ms',
           field='Vy Cma',
           spw='0,1',
           spwmap=[[0,1],[0,1]],
           gaintable=['X1de2_VYCMa_contcal_325_cont.ph2','X1de2_VYCMa_contcal_325_cont.ap2'],
           calwt = False,        # try true?  NB flagging failed
           flagbackup = False)

# 
  os.system('rm -rf X1de2_VYCMa_contcal_325_cont_ap2.clean*')
  clean(vis = 'X1de2_VYCMa_contcal_325.ms',
        imagename='X1de2_VYCMa_contcal_325_cont_ap2.clean',
        field='Vy Cma', 
        spw=contchans,
        mode='mfs',
        multiscale=[0,4,8,12] ,
        cell='0.01arcsec',npercycle=100,imsize=512,
        psfmode='hogbom',
        mask=['ellipse[[250pix,280pix],[70pix,67pix],90deg]','ellipse[[270pix,310pix],[95pix,42pix],20deg]'],
        interactive=cleaninteractive,
        niter=2000)

# Get image statistics
  imin='X1de2_VYCMa_contcal_325_cont_ap2.clean.image'
  noise=imstat(imagename=imin,
               box='250,50,500,200')['rms'][0]

  peak=imstat(imagename=imin,
              box='200,200,330,330')['max'][0]

  print 'rms %8.4f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)   
# rms 0.0009, peak  0.218, snr   243


#####################################
#  Image whole cube at lower resolution to identify line-free channels 
# This step can be omitted if you trust our selection below.
mystep = 7
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  for s in ['0','1']:
      os.system('rm -rf X1de2_VYCMa_contcal_325_spw'+s+'_5ap2.clean*')
      clean(vis = 'X1de2_VYCMa_contcal_325.ms',
            imagename='X1de2_VYCMa_contcal_325_spw'+s+'_5ap2.clean',
            field='Vy Cma', 
            spw=s,
            mode='channel',width=5,
            cell='0.02arcsec',imsize=512,usescratch=True,
            psfmode='hogbom',
            interactive=cleaninteractive,
            niter=1000)

# The following is our best guess at a selection
contchans='0:5~224;365~379;425~479;615~739;1185~1249;1460~1569;1615~1699;1820~1919,1:0~214;440~494;580~609;750~789;830~1049;1305~1364;1490~1529;1610~1659;1775~1829;1815~1834;1895~1919'


##############################################
# Image continuum  with mfs
mystep = 8
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf X1de2_VYCMa_contcal_325_contap2.clean*')
  clean(vis = 'X1de2_VYCMa_contcal_325.ms',
        imagename='X1de2_VYCMa_contcal_325_contap2.clean',
        field='Vy Cma', 
        spw=contchans,
        mode='mfs',
        multiscale=[0,4,8,12] ,
        mask=['ellipse[[250pix,280pix],[70pix,67pix],90deg]','ellipse[[270pix,310pix],[95pix,42pix],20deg]'],
        cell='0.01arcsec',imsize=512,
        psfmode='hogbom',interactive=cleaninteractive)

# OMIT
#  exportfits(imagename='X1de2_VYCMa_contcal_325_contap2.clean.image',
#             fitsimage='X1de2_VYCMa_contcal_325_contap2.clean.fits')

#####################################
#  Subtract continuum 
mystep = 9
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf X1de2_VYCMa_contcal_325.ms.contsub')
  uvcontsub(vis = 'X1de2_VYCMa_contcal_325.ms',field='V*',fitorder=1,
            combine='spw',spw='',fitspw=contchans)


  os.system('rm -rf X1de2_VYCMa_contcal_325_1_636_ap2_sub.clean*')
  clean(vis = 'X1de2_VYCMa_contcal_325.ms.contsub',
      imagename='X1de2_VYCMa_contcal_325_1_636_ap2_sub.clean',
      field='Vy Cma', 
      spw='1:636',
      antenna='!DA50',
      mode='channel',start=636,nchan=1,
      cell='0.01arcsec',npercycle=10,imsize=512,
      interactive=cleaninteractive,
      niter=1000)


#----------------------------------------------------------------------------------
#----- Image maser and line emission and make moment maps -------------------------
#----------------------------------------------------------------------------------

# Some lines from  Kaminski et al. 2013 ApJS 209 38 are imaged below.
# However, note that it is possible that the named line may be blended with other lines,
# especially for the thermal lines. Conversely, it is possible (especially
# for the maser) that faint wings exceed the range imaged.
#

#####################################
#  Image H2O  	5(1,5)-4(2,2) 	 line spw 0
mystep = 10
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]
# Line frequencies taken from Kaminski et al. 2013 ApJS 209 38

# High resolution for maser
  os.system('rm -rf X1de2_VYCMa_contcal_325.15_H2O.clean*')
  clean(vis = 'X1de2_VYCMa_contcal_325.ms.contsub',
        imagename='X1de2_VYCMa_contcal_325.15_H2O.clean',
        field='Vy Cma', 
        spw='0',
        mode='velocity',start='-38km/s',nchan=134,
        restfreq='325.152919GHz',
        cell='0.01arcsec',imsize=512,
        psfmode='hogbom',
        weighting='briggs',robust=0.5,
        threshold='10mJy',
        mask=['ellipse[[218pix,255pix],[64pix,52pix],90deg]','ellipse[[281pix,359pix],[140pix,70pix],20deg]'],
        interactive=cleaninteractive,
        usescratch=True,
        niter=1000)

# Make moment map
  thisimage='X1de2_VYCMa_contcal_325.15_H2O.clean.image'
  chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='3')
  rms1= chanstat['rms'][0]
  chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='130')
  rms2= chanstat['rms'][0]
  rms=0.5*(rms1+rms2)
  print 'rms in a channel = '+str(rms)

  os.system('rm -rf '+thisimage+'.mom0')
  immoments(imagename = thisimage,
         moments = [0],
         axis = 'spectral',
         includepix = [rms*2,1000.],
         outfile = thisimage+'.mom0')

# Export as fits OMITTED
#  exportfits(imagename=thisimage+'.mom0',
#           fitsimage=thisimage+'.mom0.fits')

#  exportfits(imagename='X1de2_VYCMa_contcal_325.15_H2O.clean.image',
#             fitsimage='X1de2_VYCMa_contcal_325.15_H2O.clean.fits')

# Note that since the dynamic range of the maser cube is 500~1000 around the peak,
# implying a noise level more than an order of magnitude
# higher than that for the fainter line wings, we don't make a first moment map.
# However, it may be possible  to obtain a meaningful first moment map with 
# additional data selection and/or blanking

#####################################
#  Image NaCl 24-23  line spw 1
mystep = 11
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]
# Line frequencies taken from Kaminski et al. 2013 ApJS 209 38

# Lower resolution for weaker lines
  os.system('rm -rf X1de2_VYCMa_contcal_312.11_NaCl.clean*')
  clean(vis = 'X1de2_VYCMa_contcal_325.ms.contsub',
        imagename='X1de2_VYCMa_contcal_312.11_NaCl.clean',
        field='Vy Cma', 
        spw='1',
        mode='velocity',
        start='-8km/s',nchan=30,width='2km/s',
        restfreq='312.1099004GHz',
        cell='0.01arcsec',imsize=512,
        psfmode='hogbom',
        threshold='4mJy',
        mask='ellipse[[284pix,263pix],[150pix,115pix],20deg]',
        interactive=cleaninteractive,
        usescratch=True,
        niter=1000)

# Make moment maps
  thisimage='X1de2_VYCMa_contcal_312.11_NaCl.clean.image'
  chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='5')
  rms1= chanstat['rms'][0]
  chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='24')
  rms2= chanstat['rms'][0]
  rms=0.5*(rms1+rms2)
  print 'rms in a channel = '+str(rms)

  os.system('rm -rf '+thisimage+'.mom0')
  immoments(imagename = thisimage,
         moments = [0],
         axis = 'spectral',
         includepix = [rms*3,100.],
         outfile = thisimage+'.mom0')

  os.system('rm -rf '+thisimage+'.mom1')
  immoments(imagename = thisimage,
         moments = [1],
         axis = 'spectral',
         includepix = [rms*4,100.],
         outfile = thisimage+'.mom1')

# Export as fits OMITTED
#  exportfits(imagename=thisimage+'.mom0',
#           fitsimage=thisimage+'.mom0.fits')

#  exportfits(imagename=thisimage+'.mom1',
#           fitsimage=thisimage+'.mom1.fits')

#  exportfits(imagename='X1de2_VYCMa_contcal_312.11_NaCl.clean.image',
#             fitsimage='X1de2_VYCMa_contcal_312.11_NaCl.clean.image.fits')

#----------------------------------------------------------------------------------
#----- End of imaging script.
#----------------------------------------------------------------------------------