# Script for imaging, self-calibration and some analysis of ALMA
# observations of TWHya.  See
# http://www.e-merlin.ac.uk/data_red/CASA/ALMA_IYAS.html for details
# and also the credits for the TW Hya data.
# amsr@jb.man.ac.uk Oct 2015

# Use the most recent released CASA version, 4.2 or later

# DO NOT LOOK AT THIS SCRIPT until you have tried to fill in the
# parameters in twhya.py yourself.

# When you get to imstat in the imaging step: suggested values

off= '170, 10, 239, 239'
on=  '100, 100, 149, 149'

# See ** below for filled-in values ** 

thesteps = []
mysteps = []
step_title = {1: 'Inspect data',
              2: 'Flag bad data',
              3: 'First clean image',
              4: 'Self-calibration: phase',
              5: 'Self-calibration: amplitude and phase',
              6: 'Split out the self-calibrated target',
              7: 'Subtract the continuum',
              8: 'Image the spectral cube',
              9: 'Image analysis' 
}
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

### 1) Inspect data
mystep = 1
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# List data
    listobs(vis='sis14_twhya_calibrated_flagged.ms')

# Plot data
    plotms(vis='sis14_twhya_calibrated_flagged.ms',
           xaxis='uvdist',
           yaxis='amp',
           avgchannel='384',
           avgtime='30s',
           field='Ceres, TW Hya, J1037-295',
           coloraxis='field',
           plotfile='Ceres_TWHya_J1037-295_uvplot.png',
           overwrite=T)

### 2) Flag data
mystep = 2
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]


# Plot just the target on short baselines
    plotms(vis='sis14_twhya_calibrated_flagged.ms',
           xaxis='time',
           field='TW Hya',
           uvrange='<50',
           coloraxis='baseline',
           plotfile='TWHya_short_uvplot.png',
           overwrite=T)

    flagdata(vis='sis14_twhya_calibrated_flagged.ms', 
             antenna='DV19',
             scan='16')

# enter the other data to be flagged
    flagdata(vis='sis14_twhya_calibrated_flagged.ms', 
             antenna='DV20',  #**
             scan='36')       #**

### 3) First clean image
mystep = 3
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    clearcal('sis14_twhya_calibrated_flagged.ms') # Remove previous self-cal
    os.system('rm -rf twhya_cont_clean0*') #   Delete any old images of the same name
    clean(vis='sis14_twhya_calibrated_flagged.ms',
          imagename='twhya_cont_clean0',
          field='TW Hya',
          mode='mfs',
          imsize=[250,250],          # ** image size in pixels
          cell=['0.08arcsec'],       # ** pixel size
          weighting='natural',
          interactive=True)

# ** Make sure you have entered xoff etc. near start of script **

    rms1=imstat(imagename='twhya_cont_clean0.image',
                box=off)['rms'][0]  

    peak1=imstat(imagename='twhya_cont_clean0.image',
                 box=on)['max'][0]     

    print 'twhya_cont_clean0.image rms %7.3f mJy' % (rms1*1000.)
    print 'twhya_cont_clean0.image peak %7.3f mJy/bm' % (peak1*1000.)
    print 'twhya_cont_clean0.image S/N %7.f' % (peak1/rms1)

### 4) Self-calibration: phase
mystep = 4
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# Plot antenna locations to select refant
    plotants(vis='sis14_twhya_calibrated_flagged.ms',
         figfile='sis14_twhya_ants.png')




# Plot the model (FT of the CC from the image you just made)
    plotms(vis='sis14_twhya_calibrated_flagged.ms',
           xaxis='uvdist', 
           yaxis='amp', 
           ydatacolumn='model',
           avgchannel='384',
           avgtime='30s',
           field='TW Hya',
           plotfile='TWHya_model_uvplot.png',
           overwrite=T)

# The data themselves
    plotms(vis='sis14_twhya_calibrated_flagged.ms',
           xaxis='uvdist', 
           yaxis='amp', 
           ydatacolumn='data',
           avgchannel='384',
           avgtime='30s',
           field='TW Hya',
           coloraxis='baseline',
           plotfile='TWHya_data_uvplot.png',
           overwrite=T)

# Phase calibration
    os.system('rm -rf twhya_cont.ph1')
    gaincal(vis='sis14_twhya_calibrated_flagged.ms',
            caltable='twhya_cont.ph1',
            field='TW Hya',
            solint='30s',
            calmode='p',        #**
            refant='DV22')      #**

# Plot the solutions
    plotcal(caltable='twhya_cont.ph1',  #**
            xaxis='time',                       #**
            yaxis='phase',                      #**
            subplot=531,
            iteration='antenna',
            plotrange=[-1,-1,-180,180],
            figfile='sis14_p1_phase_scan.png')

# Apply the solutions
    applycal(vis='sis14_twhya_calibrated_flagged.ms',
             field='TW Hya',
             gaintable=['twhya_cont.ph1'])

    os.system('rm -rf twhya_cont_clean_p1*') #   Delete any old images of the same name
    clean(vis='sis14_twhya_calibrated_flagged.ms',
          imagename='twhya_cont_clean_p1',
          field='TW Hya',                  #**
          mode='mfs',                      #**
          imsize=[250,250],                #**
          cell=['0.08arcsec'],             #**
          weighting='natural',             #**
          niter=1000,
          npercycle=200,
          interactive=True)

# ** Make sure you have entered xoff etc. near start of script **

    rms1=imstat(imagename='twhya_cont_clean_p1.image',
                box=off)['rms'][0]  
    
    peak1=imstat(imagename='twhya_cont_clean_p1.image',
                 box=on)['max'][0]     


    print 'twhya_cont_clean_p1.image rms %7.3f mJy' % (rms1*1000.)
    print 'twhya_cont_clean_p1.image peak %7.3f mJy/bm' % (peak1*1000.)
    print 'twhya_cont_clean_p1.image S/N %7.f' % (peak1/rms1)

### 5) Self-calibration: amplitude and phase 
mystep = 5
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf twhya_cont.ap1')
    gaincal(vis='sis14_twhya_calibrated_flagged.ms',
            caltable='twhya_cont.ap1',
            field='TW Hya',
            solint='30s',
            calmode='ap',
            gaintable='twhya_cont.ph1',
            refant='DV22')

    plotcal(caltable='twhya_cont.ap1',
            xaxis='time',
            yaxis='phase',
            subplot=531,
            iteration='antenna',
            plotrange=[-1,-1,-180,180],
            figfile='sis14_ap1_phase_scan.png')
    
    plotcal(caltable='twhya_cont.ap1',
            xaxis='time',
            yaxis='amp',
            subplot=531,
            iteration='antenna',
            plotrange=[-1,-1,0,2],
            figfile='sis14_ap1_amp_scan.png')

# Apply the solutions
    applycal(vis='sis14_twhya_calibrated_flagged.ms',
             field='TW Hya',
             gaintable=['twhya_cont.ph1','twhya_cont.ap1'])

# Image again
    os.system('rm -rf twhya_cont_clean_ap1*') #   Delete any old images of the same name
    clean(vis='sis14_twhya_calibrated_flagged.ms',
          imagename='twhya_cont_clean_ap1',
          field='TW Hya',               #**
          mode='mfs',                   #**
          imsize=[250,250],             #**
          cell=['0.08arcsec'],          #**
          weighting='natural',          #**
          niter=1000,                   #**
          npercycle=200,                #**
          interactive=True)




    rms1=imstat(imagename='twhya_cont_clean_ap1.image',
                box=off)['rms'][0]  

    peak1=imstat(imagename='twhya_cont_clean_ap1.image',
                 box=on)['max'][0]     

    print 'twhya_cont_clean_ap1.image rms %7.3f mJy' % (rms1*1000.)
    print 'twhya_cont_clean_ap1.image peak %7.3f mJy/bm' % (peak1*1000.)
    print 'twhya_cont_clean_ap1.image S/N %7.f' % (peak1/rms1)


### 6) Split out the self-calibrated target
mystep = 6
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf twhya_calibrated_flagged.ms')
    split(vis='sis14_twhya_calibrated_flagged.ms',
          outputvis='twhya_calibrated_flagged.ms',
          field='TW Hya',
          keepflags=F)


### 7) Subtract the continuum
mystep = 7
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    plotms(vis='twhya_calibrated_flagged.ms',
           xaxis='channel',
           yaxis='amp',
           ydatacolumn='data',
           uvrange='0~50',
           avgtime='1e9',
           avgscan=True,
           avgbaseline=True,
           plotfile='TWHya_cal_spectrum.png',
           overwrite=T)

    os.system('rm -rf twhya_calibrated_flagged.ms.contsub')
    uvcontsub(vis = 'twhya_calibrated_flagged.ms',
              field = 'TW Hya',
              fitspw = '0:240~280',   # **
              excludechans = True,
              fitorder = 1,
              solint='int')

### 8) Image the spectral cube
mystep = 8
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf twhya_n2hp.*')
    clean(vis = 'twhya_calibrated_flagged.ms.contsub',
          imagename = 'twhya_n2hp',
          field = 'TW Hya',
          mode = 'velocity',
          nchan = 15,
          start = '0.0km/s',
          width = '0.5km/s',
          outframe = 'LSRK',
          restfreq = '372.67249GHz',
          interactive = T,
          imsize = [250, 250],
          cell = '0.08arcsec',
          weighting = 'briggs',
          robust = 0.5)


### 9) Image analysis
mystep = 9
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf twhya_n2hp.mom0')
    immoments(imagename='twhya_n2hp.image',
              outfile='twhya_n2hp.mom0',
              includepix=[20e-3,100],      #**
              chans='4~12',                #**
              moments=[0])


    imview(raster={'file': 'twhya_n2hp.mom0',
                   'colorwedge': T},
           contour={'file': 'twhya_cont_clean_ap1.image',  
                    'unit': 0.004,
                    'levels': [-1,1,4,16,64] },  
           zoom=2,
           out='twhya_n2hp.mom0_cont.png')  
    

    exportfits(imagename='twhya_n2hp.image',
               fitsimage='twhya_n2hp.image.fits')