# Script to image eMERLIN data with CASA 4.2.2 or later

# 1252+5634.ms contains target 3C277.1 with all external calibration
# applied, including the phase-reference solutions

target = '1252+5634' # 3C277.1
antref='Mk2'          # refant - usually - change if known to be bad.

# Enter step(s) as mysteps

# Fill in the missing parameters and values

#############################
thesteps = []
step_title = {1: 'Check phase-reference calibrated target data',
              2: 'Make first target image',
              3: 'Choose interval for first phase self-cal',
              4: 'First phase self-cal',
              5: 'Apply solutions and image',
              6: 'Second phase self-cal',
              7: 'Apply new phase calibration and image',
              8: 'Chose interval for amp self-cal',
              9: 'Amplitude self-calibration',
              10: 'Apply amp self-cal and image',
              11: 'Reweight according to sensitivity',
              12: 'Multiscale clean',
              13: 'Clean with uvtaper',
              14: 'Clean with Briggs robust weighting'}

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) Check target MS
mystep = 1

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf '+target+'.ms.listobs')
    listobs(vis=target+'.ms',
            listfile=target+'.ms.listobs')

# plot amp against uv distance            
    plotms(vis=target+'.ms',
           xaxis='**',
           yaxis='**',
           correlation='**',
           avgchannel='**',
           coloraxis='spw',
           plotfile='3C277.1_uvdist.png', 
           overwrite=True)


### 2) Image target
# Use known properties of 3C277.1 as guidance in initial masking.
# Expand mask as more emission appears above decreasing noise 
# in successive cycles.

mystep = 2

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]


    os.system('rm -rf '+target+'_1.clean*')
    clean(vis=target+'.ms',
          imagename=target+'_1.clean',
          cell='**arcsec',
          niter=700,
          imsize='***',
          interactive=True)

    rms1=imstat(imagename=target+'_1.clean.image',
                  box='10,400,199,499')['rms'][0]
    peak1=imstat(imagename=target+'_1.clean.image',
                  box='10,10,499,499')['max'][0]
 
    print target+'_1.clean.image rms %7.3f mJy' % (rms1*1000.)
    print target+'_1.clean.image peak %7.3f mJy/bm' % (peak1*1000.)


    imview(raster={'file': target+'_1.clean.image',
                   'colormap':'Rainbow 2',
                   'scaling': -2,
                   'range': [-1.*rms1, 10.*rms1],
                   'colorwedge': True},
           zoom=2,
           out=target+'_1.clean.png')

    plotms(vis=target+'.ms',
           xaxis='uvdist',
           yaxis='amp',
           ydatacolumn='model',
           correlation='RR,LL',
           avgchannel='64',
           coloraxis='spw',
           plotfile='3C277.1_uvdist_model1.png', 
           overwrite=True)


### 3) Chose interval for first phase self-cal
# If you like, inspect  phases against time to help decide solint.
# However, phase slopes are also due to source structure.
# Look at the shortest baseline (Mk2&Pi) 
# But, inadvisable to use a very short solint while model imperfect.
# Half a scan length seems OK.

mystep = 3

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    plotms(vis=target+'.ms',
           xaxis='time',
           yaxis='phase',
           antenna='Mk2&Pi',
           correlation='RR,LL',
           avgchannel='64',
           coloraxis='spw',
           plotfile='3C277.1_phase_1.png', 
           overwrite=True)

### 4) First phase self-cal

mystep = 4

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf target.p1')
    gaincal(vis=target+'.ms',
            calmode='***',
            caltable='target.p1',
            field=target,
            solint='***',
            refant=antref,
            minblperant=2,minsnr=2)
 
# Plot solutions
    plotcal(caltable='target.p1',xaxis='time',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        poln='L',figfile='target.p1_L.png')

    plotcal(caltable='target.p1',xaxis='time',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        poln='R',figfile='target.p1_R.png')

### 5) Apply phase solutions and image again

mystep = 5

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    applycal(vis=target+'.ms',
             gaintable='***')


# clean phase self-calibrated image
    os.system('rm -rf '+target+'_p1.clean*')
    clean(vis=target+'.ms',
          imagename=target+'_p1.clean',
          cell='***arcsec',
          niter=1000, # stop before this if residuals reach noise
          imsize='***',
          interactive=True)

    rmsp1=imstat(imagename=target+'_p1.clean.image',
                  box='10,400,199,499')['rms'][0]
    peakp1=imstat(imagename=target+'_p1.clean.image',
                  box='10,10,499,499')['max'][0]
 
    print target+'_p1.clean.image rms %7.3f mJy' % (rmsp1*1000.)
    print target+'_p1.clean.image peak %7.3f mJy/bm' % (peakp1*1000.)


    viewer(target+'_p1.clean.image')

    imview(raster={'file': target+'_p1.clean.image',
                   'range': [-1.*rmsp1, 20.*rmsp1],
                   'colorwedge': True},
           zoom=2,
           out=target+'_p1.clean.png')

# Has SNR improved?

### 6) Self-cal target #2
# Better model, shorter solution interval
# This time, make spectral index image to prepare for amp self-cal

mystep = 6

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf target.p2')
    gaincal(vis=target+'.ms',
            calmode='***',
            caltable='target.p2',
            field=target,
            solint='***',
            refant=antref,
            minblperant=2,minsnr=2)
 
    # Plot solutions
    plotcal(caltable='target.p2',xaxis='time',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        poln='L',figfile='target.p2_L.png')

    plotcal(caltable='target.p2',xaxis='time',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        poln='R',figfile='target.p2_R.png')


### 7) Apply new phase calibration and image
# This time, make spectral index image to prepare for amp self-cal

mystep = 7

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    applycal(vis=target+'.ms',
             gaintable='***')

# clean phase self-calibrated image
    os.system('rm -rf '+target+'_p2.clean*')
    clean(vis=target+'.ms',
          imagename=target+'_p2.clean',
          nterms=2, # Make spectral index image
          cell='***arcsec',
          niter=1000, # stop before this if residuals reach noise
          imsize='***',
          interactive=True)

    rmsp2=imstat(imagename=target+'_p2.clean.image.tt0',
                  box='10,400,199,499')['rms'][0]
    peakp2=imstat(imagename=target+'_p2.clean.image.tt0',
                  box='10,10,499,499')['max'][0]
 
    print target+'_p2.clean.image.tt0 rms %7.3f mJy' % (rmsp2*1000.)
    print target+'_p2.clean.image.tt0 peak %7.3f mJy/bm' % (peakp2*1000.)


    viewer(target+'_p2.clean.image.tt0')

    imview(raster={'file': target+'_p2.clean.image.tt0',
                   'range': [-1.*rmsp2, 20.*rmsp2],
                   'colorwedge': True},
           zoom=2,
           out=target+'_p2.clean.png')


### 8) Chose interval for amplitude self-cal

mystep = 8

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]


    plotms(vis=target+'.ms',
           xaxis='time',
           yaxis='amp',
           ydatacolumn='corrected',
           antenna='Mk2&Pi',
           correlation='RR,LL',
           avgchannel='64',
           coloraxis='spw',
           plotfile='3C277.1_amp_1.png', 
           overwrite=True)


### 9) amplitude self-cal

mystep = 9

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf target.ap1')
    gaincal(vis=target+'.ms',
            calmode='***',              # amp and phase
            caltable='target.ap1',
            field=target,
            solint='***',
            refant=antref,
            gaintable='***',
            minblperant=2,minsnr=2)
 
    # Plot solutions
    # If these are all more than ~20% from unity, check for a poor model
    # If just a few are  more than ~20% from unity, check for bad data
    plotcal(caltable='target.ap1',xaxis='time',
        iteration='antenna',yaxis='amp',subplot=321,
        poln='L',figfile='target.ap1_L.png')

    plotcal(caltable='target.ap1',xaxis='time',
        iteration='antenna',yaxis='amp',subplot=321,
        poln='R',figfile='target.ap1_R.png')


### 10) Apply amplitude calibration and image

mystep = 10

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    applycal(vis=target+'.ms',
             gaintable=['target.p2','target.ap1'])

# clean amp&phase self-calibrated image
    os.system('rm -rf '+target+'_ap1.clean*')
    clean(vis=target+'.ms',
          imagename=target+'_ap1.clean',
          nterms='**',
          cell='***arcsec',
          niter=1500, # stop before this if residuals reach noise
          imsize='***',
          interactive=True)

    rmsap1=imstat(imagename=target+'_ap1.clean.image.tt0',
                  box='10,400,199,499')['rms'][0]
    peakap1=imstat(imagename=target+'_ap1.clean.image.tt0',
                  box='10,10,499,499')['max'][0]
 
    print target+'_ap1.clean.image.tt0 rms %7.3f mJy' % (rmsap1*1000.)
    print target+'_ap1.clean.image.tt0 peak %7.3f mJy/bm' % (peakap1*1000.)


    viewer(target+'_ap1.clean.image.tt0')

    imview(raster={'file': target+'_ap1.clean.image.tt0',
                   'range': [-1.*rmsap1, 20.*rmsap1],
                   'colorwedge': True},
           zoom=2,
           out=target+'_ap1.clean.png')

# 
# 1252+5634_ap1.clean.image.tt0 rms   0.107 mJy
# 1252+5634_ap1.clean.image.tt0 peak 172.318 mJy/bm


### 11) Re-weight according to sensitivity
# Re-weight the data according to sensitivity (rms in short intervals)

mystep = 11

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# View weights
    plotms(vis=target+'.ms', field=target, xaxis='time', 
           yaxis='weight',antenna=antref+'&*',
           correlation='LL,RR',
           showgui=gui,
           coloraxis='antenna2',plotfile='target_weight-time.png',
           overwrite=True)

# The weights may look odd because e-MERLIN weights are not handled
# correctly by all versions of CASA

# re-weight
    statwt(vis=target+'.ms')

# View weights again
    plotms(vis=target+'.ms', field=target, xaxis='uvdist', 
           yaxis='weight',antenna=antref+'&*',
           correlation='LL,RR',
           showgui=gui,
           coloraxis='antenna2',
           plotfile='target_re-weight-uv.png',
           overwrite=True)


### 12) Multiscale clean
# Try searching for multiple scales in clean

mystep = 12

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# clean amp&phase self-calibrated image
    os.system('rm -rf '+target+'_ap1_multi.clean*')
    clean(vis=target+'.ms',
          imagename=target+'_ap1_multi.clean',
          nterms='***',
          cell='***arcsec',
          niter='***', # 
          multiscale=['**','**','**'],
          imsize='***',
          interactive=True)
# More components are subtracted per cycle so clean runs faster


    rmsmulti=imstat(imagename=target+'_ap1_multi.clean.image.tt0',
                  box='10,400,199,499')['rms'][0]
    peakmulti=imstat(imagename=target+'_ap1_multi.clean.image.tt0',
                  box='10,10,499,499')['max'][0]
 
    print target+'__ap1_multi.clean.image.tt0 rms %7.3f mJy' % (rmsmulti*1000.)
    print target+'__ap1_multi.clean.image.tt0 peak %7.3f mJy/bm' % (peakmulti*1000.)


### 13) Different imaging weights: tapering

mystep = 13

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# Experiment with changing the relative weights given to longer
# baselines.  Applying a taper reduces this, giving a larger
# synthesised beam more sensitive to extended emission. This just
# changes the weights during gridding for imaging, no permanent change
# to the uv data. However, strong tapering increases the noise.

# Look at the uv distance plot and choose an outertaper (half-power of
# Gaussian taper applied to data) which is about 2/3 of the longest
# baseline. Units are wavelengths, so calculate that or plot amp
# v. uvwave

    os.system('rm -rf '+target+'_ap1_taper.clean*')
    clean(vis=target+'.ms',
          imagename=target+'_ap1_taper.clean',
          cell='***arcsec',
          nterms='***',
          niter=1500,
          multiscale=['**','**','**'],
          uvtaper=True,
          outertaper='***',
          imsize='***',
          interactive=True)

    rmstaper=imstat(imagename=target+'_ap1_taper.clean.image.tt0',
                  box='10,400,199,499')['rms'][0]
    peaktaper=imstat(imagename=target+'_ap1_taper.clean.image.tt0',
                  box='10,10,499,499')['max'][0]
 
    print target+'_ap1_taper.clean.image.tt0 rms %7.3f mJy' % (rmstaper*1000.)
    print target+'_ap1_taper.clean.image.tt0 peak %7.3f mJy/bm' % (peaktaper*1000.)


### 14) Different imaging weights: Briggs robust

mystep = 14

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# Experiment with changing the relative weights given to longer
# baselines.  Default natural weighting gives zero weight to parts of
# the uv plane with no samples.  Briggs weighting interpolates between
# samples, giving more weight to the interpolated areas for lower
# values of robust.  For an array like e-MERLIN with sparser sampling
# on long baselines, this has the effect of increasing their
# contribution, i.e. making the synthesised beam smaller. This may
# bring up the jet/core, at the expense of the extended lobes.  Low
# values of robust increase the noise. For now, remove tapering and
# multiscale.

    os.system('rm -rf '+target+'_ap1_robust.clean*')
    clean(vis=target+'.ms',
          imagename=target+'_ap1_robust.clean',
          cell='***arcsec',
          nterms='***',
          niter='***',
          weighting='***',
          robust='***',
          imsize='***',
          interactive=True)

    rmsrobust=imstat(imagename=target+'_ap1_robust.clean.image.tt0',
                  box='10,400,199,499')['rms'][0]
    peakrobust=imstat(imagename=target+'_ap1_robust.clean.image.tt0',
                  box='10,10,499,499')['max'][0]
 
    print target+'_ap1_robust.clean.image.tt0 rms %7.3f mJy' % (rmsrobust*1000.)
    print target+'_ap1_robust.clean.image.tt0 peak %7.3f mJy/bm' % (peakrobust*1000.)