## Reducing HD data

#
# step 0: Get the visibility data and run listobs and plotants
# step 1: Put in flux density of phase reference source and store flags
# step 2: Additional flags (which come from the next step, so don't
#         put them in at first.
#         Then run gaincal and plot ampl and phases.  This will suggest
#         the additional flags, add in then run again.
# Step 3: Apply flags again (if beginning here)
#         Apply the calibration from the phase reference source
# Step 4: tclean using normal single resolution clean.
#         Do interactively.  Can try usemask = 'auto-multithresh' is you want.
# Step 4.1: Does contour/raster plot
#
# Step 6: Do phase selfcal using 60-s solution interval
# Step 6.1: Does phase solution plots.  Noisish, but okay
# Step 7: Apply this self-cal  7.1 plots
# Step 8: tclean the self-cal data, single resolution  (8.1 to plot)
# Step 9: tclean the self-cal data using nominal multi-scale
# Step 10: tclean the self-cal data using point source weight multi-scale
# Step 16,16.1,17,17.1  gain selfcal
#         
zvis = 'HD.ms'
#
dostep = [8.1]
#
#
#
if 0 in dostep:
    # Copy the data from disk.
    listobs(vis=zvis,listfile=zvis+'.listobs')
    # contains 2 SB's.  All calibrations done except for phase referencing
#
if 1 in dostep:
    print 'put in flux density of phase-ref cal, and store original flags'
    setjy(vis=zvis,field='J1742-1517',usescratch=True,standard='manual')
    flagmanager(vis=zvis, mode='save',versionname='Original')
#
if 2 in dostep:
    print 'flag and do gaincal of phase cal'
    flagdata(vis=zvis,spw = '3', flagbackup = False)
    flagdata(vis=zvis,antenna='DV03', flagbackup = False)
    os.system('rm -rf pref')
    gaincal(vis=zvis,
            field ='J1742-1517',
            caltable = 'pref',
            spw = '0,1,2',
            solint = 'inf',
            refant = 'DA51')
    plotcal(caltable = 'pref',
            xaxis = 'scan',
            yaxis = 'amp',
            subplot = 311,
            plotrange = [0,0,0,0.01],
            iteration = 'spw,antenna')
    plotcal(caltable = 'pref',
            xaxis = 'scan',
            yaxis = 'phase',
            spw = '0', poln = 'X',
            subplot = 311,
            iteration = 'antenna',
            plotrange=[0,0,-180,180])
    # DV03 baseline error
#
if 3 in dostep:
    print 'apply to target with flagging'
    flagmanager(vis=zvis,mode='restore',versionname='Original')
    flagdata(vis=zvis,spw = '3', flagbackup = False)
    flagdata(vis=zvis,antenna='DV03,DA65', flagbackup = False)
    applycal(vis=zvis,
             gaintable = 'pref',
             spw = '0,1,2',
             flagbackup = False)
#
if 4 in dostep:
    print 'image target'
    Im = 'HD_2F'
    os.system('rm -rf '+Im+'*')
    tclean(vis=zvis,
           field = 'HD163296',
           intent = '*OBSERVE_TAR*',
           imagename=Im, 
           cell = '0.04arcsec', 
           imsize = 500,
           weighting = 'briggs',
           robust = 0.5,
           restart = True,
           mask = 'HD_2.mask',
           interactive = True,
           deconvolver = 'clark',
           niter = 10000,
           usemask = 'user')
           #mask = 'HD_2.mask')
           #usemask = 'auto-thresh2',
           #maskthreshold = 3.0,
           #maskresolution = '0,2arcsec',
           #nmask = 5,
           #pbmask = 0.2)
           #usemask = 'auto-multithresh',
           #noisethreshold = 3.0,
           #sidelobethreshold = 3.0)
#
if 4.1 in dostep:
    Im = 'HD_2F'
    Jm = Im+'.image'
    imax = imstat(Jm)['max'][0]
    imin = imax/50.0
    imview (raster={'file': Jm,
                    'range': [-imin, imax],
                    'colormap': 'Smooth 1', 'colorwedge': True},
            contour={'file': Jm,
                     'levels': [-1,1,2,3,4,8,16,24,25,26,35,45],
                     'unit': float(imin)},
            zoom={'blc':[186,186],'trc':[326,326]}) 

if 6 in dostep:
    print 'do the phase selfcal'
    # pscalF  both pols, 60 sec.
    print 'redo flagging'
    flagmanager(vis=zvis,mode='restore',versionname='Original')
    flagdata(vis=zvis,spw = '3', flagbackup = False)
    flagdata(vis=zvis,antenna='DV03,DA65', flagbackup = False)
    pmod = 'pscalF'
    os.system('rm -rf '+pmod)
    print 'redo making first model'
    ft(vis=zvis,field='HD163296',model='HD_2F.model')
    gaincal(vis=zvis,
            field ='HD163296',
            caltable = pmod,
            gaintable = 'pref',
            gaintype = 'T',
            calmode = 'p',
            spw = '0,1,2',
            combine = 'spw,scan',
            solint = '60s',
            minsnr = 1.0,
            refant = 'DA51')
#
if 6.1 in dostep:
    print 'plot phase selfcal'
    pmod = 'pscalF'
    plotcal(caltable = pmod,
            xaxis = 'scan',
            yaxis = 'phase',
            spw = '0', poln = 'X',
            subplot = 221,
            iteration = 'antenna',
            plotrange=[0,0,-180,180])
#
if 7 in dostep:
    print 'apply selfcal'
    applycal(vis=zvis,
             gaintable = ['pref',pmod],
             spw = '0,1,2',
             spwmap = [[0,1,2],[0,0,0]],
             flagbackup = False)
    
if 8 in dostep:
    print 'image target'
    Im = 'HD_2F'+pmod
    os.system('rm -rf '+Im+'*')
    tclean(vis=zvis,
           field = 'HD163296',
           intent = '*OBSERVE_TAR*',
           imagename=Im, 
           cell = '0.04arcsec', 
           imsize = 500,
           weighting = 'briggs',
           robust = 0.5,
           restart = True,
           interactive = True,
           deconvolver = 'clark',
           niter = 10000,
           mask = 'HD_2F.mask',
           usemask = 'user')
           #mask = 'HD_1.mask')
           #usemask = 'auto-thresh2',
           #maskthreshold = 3.0,
           #maskresolution = '0,2arcsec',
           #nmask = 5,
           #pbmask = 0.2)
           #usemask = 'auto-multithresh',
           #noisethreshold = 3.0,
           #sidelobethreshold = 3.0)
#
if 8.1 in dostep:
    Im = 'HD_2F'+pmod
    Jm = Im+'.image'
    imax = imstat(Jm)['max'][0]
    imin = imax/400.0
    imview (raster={'file': Jm,
                    'range': [-imin, imax],
                    'colormap': 'Smooth 1', 'colorwedge': True},
            contour={'file': Jm,
                     'levels': [-1,1,2,4,8,16,50,100,198,200,202,280,350],
                     'unit': float(imin)},
            zoom={'blc':[186,186],'trc':[326,326]}) 
if 9 in dostep:
    print 'image target'
    Im = 'HD_2FMS'+pmod
    os.system('rm -rf '+Im+'*')
    tclean(vis=zvis,
           field = 'HD163296',
           intent = '*OBSERVE_TAR*',
           imagename=Im, 
           cell = '0.04arcsec', 
           imsize = 500,
           weighting = 'briggs',
           robust = 0.5,
           restart = True,
           interactive = True,
           deconvolver = 'multiscale',
           scales = [0,7,15],
           niter = 10000,
           mask = 'HD_2F.mask',
           usemask = 'user')
           #mask = 'HD_1.mask')
           #usemask = 'auto-thresh2',
           #maskthreshold = 3.0,
           #maskresolution = '0,2arcsec',
           #nmask = 5,
           #pbmask = 0.2)
           #usemask = 'auto-multithresh',
           #noisethreshold = 3.0,
           #sidelobethreshold = 3.0)
#
if 9.1 in dostep:
    Im = 'HD_2FMS'+pmod
    Jm = Im+'.image'
    imax = imstat(Jm)['max'][0]
    imin = imax/400.0
    imview (raster={'file': Jm,
                    'range': [-imin, imax],
                    'colormap': 'Smooth 1', 'colorwedge': True},
            contour={'file': Jm,
                     'levels': [-1,1,2,4,8,16,50,100,198,200,202,280,350],
                     'unit': float(imin)},
            zoom={'blc':[186,186],'trc':[326,326]}) 
#
if 10 in dostep:
    print 'image target'
    Im = 'HD_2FMS3'+pmod
    os.system('rm -rf '+Im+'*')
    tclean(vis=zvis,
           field = 'HD163296',
           intent = '*OBSERVE_TAR*',
           imagename=Im, 
           cell = '0.04arcsec', 
           imsize = 500,
           weighting = 'briggs',
           robust = 0.5,
           restart = True,
           interactive = True,
           deconvolver = 'multiscale',
           scales = [0,7,15],
           smallscalebias = 0.8,
           niter = 10000,
           mask = 'HD_2F.mask',
           usemask = 'user')
           #mask = 'HD_1.mask')
           #usemask = 'auto-thresh2',
           #maskthreshold = 3.0,
           #maskresolution = '0,2arcsec',
           #nmask = 5,
           #pbmask = 0.2)
           #usemask = 'auto-multithresh',
           #noisethreshold = 3.0,
           #sidelobethreshold = 3.0)
#
if 10.1 in dostep:
    Im = 'HD_2FMS3'+pmod
    Jm = Im+'.image'
    imax = imstat(Jm)['max'][0]
    imin = imax/400.0
    imview (raster={'file': Jm,
                    'range': [-imin, imax],
                    'colormap': 'Smooth 1', 'colorwedge': True},
            contour={'file': Jm,
                     'levels': [-1,1,2,4,8,16,50,100,198,200,202,280,350],
                     'unit': float(imin)},
            zoom={'blc':[186,186],'trc':[326,326]})
#
if 16 in dostep:
    print 'do the phase and amp selfcal'
    # pscalG  both pols, 500 sec.
    pmod = 'pscalG'
    print 'put in self-cal model'
    ft(vis=zvis,model='HD_2FpscalF.model',field='HD163296')
    os.system('rm -rf '+pmod)
    gaincal(vis=zvis,
            field ='HD163296',
            caltable = pmod,
            gaintable = ['pref','pscalF'],
            spwmap = [[0,1,2],[0,0,0]],
            calmode = 'ap',
            spw = '0,1,2',
            combine = 'scan',
            solint = '500s',
            minsnr = 1.0,
            refant = 'DA51')
#
if 16.1 in dostep:
    print 'plot phase selfcal'
    pmod = 'pscalG'
    plotcal(caltable = pmod,
            xaxis = 'scan',
            yaxis = 'amp',
            #spw = '0', poln = 'X',
            subplot = 221,
            iteration = 'antenna',
            plotrange=[0,0,0.8,1.2])
#
if 17 in dostep:
    print 'apply selfcal'
    applycal(vis=zvis,
             gaintable = ['pref','pscalF','pscalG'],
             spw = '0,1,2',
             spwmap = [[0,1,2],[0,0,0],[0,1,2]],
             flagbackup = False)
    
if 18 in dostep:
    print 'image target'
    Im = 'HD_2F'+pmod
    os.system('rm -rf '+Im+'*')
    tclean(vis=zvis,
           field = 'HD163296',
           intent = '*OBSERVE_TAR*',
           imagename=Im, 
           cell = '0.04arcsec', 
           imsize = 500,
           weighting = 'briggs',
           robust = 0.5,
           restart = True,
           interactive = True,
           deconvolver = 'clark',
           niter = 10000,
           mask = 'HD_2FpscalF.mask',
           usemask = 'user')
           #mask = 'HD_1.mask')
           #usemask = 'auto-thresh2',
           #maskthreshold = 3.0,
           #maskresolution = '0,2arcsec',
           #nmask = 5,
           #pbmask = 0.2)
           #usemask = 'auto-multithresh',
           #noisethreshold = 3.0,
           #sidelobethreshold = 3.0)
#
if 18.1 in dostep:
    Im = 'HD_2F'+pmod
    Jm = Im+'.image'
    imax = imstat(Jm)['max'][0]
    imin = imax/400.0
    imview (raster={'file': Jm,
                    'range': [-imin, imax],
                    'colormap': 'Smooth 1', 'colorwedge': True},
            contour={'file': Jm,
                     'levels': [-1,1,2,4,8,16,50,100,198,200,202,280,350],
                     'unit': float(imin)},
            zoom={'blc':[186,186],'trc':[326,326]})