# MWA School raw script for you to correct and improve

# Analysis of NGC 4258 = 12190+47182 
# K band continuum analysis

# This script is supposed to analyse a part of the VLA dataset AH847_1 which you have already seen
# in the previous example ngc4258-lines.py.
# However, it is not finished! Please fix the following problems:
# (you can look for them problems by finding "_here" with your editor)

# 1) Determine the spectral window IDs of the continuum and enter them in the split call to
#    create a separate MS with only the data interesting for this analysis
#
# 2) Learn about the flagdata task and understand the "quacking" of VLA data.
#    Determine the correct quack interval and enter the number in the corresponding flagdata task.
#
# 3) Check the MS with viewer: see that in LL, LR, and RL, there was a problem in the middle of the time range
#    lasting 9 integrations start is 2:48:25. Determine the end time and enter it in the flagdata task.
#    (Hint: use the Scaling Power Cycles under Basic Settings in the Data Display Options panel
#     to adjust the intensity scaling, also changing the color map to Hot Metal 1 helps.
#     See the exact time by zooming in and moving the pointer over the cell in question. Then look
#     at the text field at the bottom of the Viewer Display Panel.) 
#
# 4) Learn about the setjy task and find to what value you have to set the "fluxdensity" parameter 
#    in order to determine the fluxdensity value from the internal reference table.
#
# 5) Learn about the fluxscale task and find out what the reference parameter means and how you
#    have to set it in the context of this analysis to transfer the flux density scale from the
#    flux calibrator to the phase calibrator
#
# 6) Learn about the applycal task and its gaincurve parameter. Find what option it should be set to
#    for this dataset.
#
# 7) Learn about the clean task and complete the various missing parameters in the clean call
#    for cleaning the phase calibrator
#
# 8) Look at the image of the phase calibrator and determine the size of the synthesized beam using
#    the imhead task, mode = list.
#    (Note: In the logger output of imhead, beammajor and beamminor give the major and minor
#    diameter of the beam ellipse in arcsecs.)
#
# 9) Choose the restoring beam in the clean call for the target to be a circular beam with the same area.
#    Choose the cell parameter about 3 to 5 times smaller.
#
# 10)Now add a call to the imfit task to fit a gaussian in the box '54, 54, 74, 74'
#    to the target image 


doall = False
try:
    print "List of steps to be executed ... ", mysteps
except:
    print "global variable mysteps not set."
    mysteps = []
if (mysteps==[]):
    print "mysteps empty. Executing all steps."
    doall = True

# my globals
basename = 'ah847_1-k-continuum'    

steptitle = {
    0 : "import",
    1 : "split out the data of interest and record original flags",
    2 : "reset flags and apply my flagging",
    3 : "clear, recalculate, and apply calibration",
    4 : "clean phase calibrator",
    5 : "determine synthesized beam size",
    6 : "clean the target",
    7 : "fit a gaussian to the image"
    }

# import
mystep = 0
if(mystep in mysteps or doall):
    print "Step ", mystep, steptitle[mystep]
    
    os.system('rm -rf '+basename+'*')

    importvla(
        archivefiles = 'AH847_1',
        vis = basename+'.ms',
        bandname = 'K',
        project = '',
        starttime = '',
        stoptime = '',
        applytsys = True,
        autocorr = False,
        antnamescheme = 'new', # better to avoid mix-up of antenna name and ID 
        keepblanks = False,
        evlabands = False,
        async = False
        )


mystep = 1
if(mystep in mysteps or doall):
    print "Step ", mystep, steptitle[mystep]
    
    os.system('rm -rf '+basename+'-selected*')

    split(
        vis = basename+'.ms',
        outputvis = basename+'-selected.ms',
        datacolumn = 'data',
        spw = mystring_here, # problem 1
        field = '1,2,4' 
        )

    # record flags
    flagmanager(
        vis = basename+'-selected.ms',
        mode = 'save',
        versionname = 'original_data',
        merge = 'replace'
        )

# flagging
mystep = 2
if(mystep in mysteps or doall):
    print "Step ", mystep, steptitle[mystep]

    # restore original flags
    flagmanager(
        vis = basename+'-selected.ms',
        mode = 'restore',
        versionname = 'original_data',
        merge = 'replace'
        )
    
    flagdata(
        vis = basename+'-selected.ms',
        mode = 'quack',
        unflag = False,
        quackinterval = mynumber_here, # problem 2 
        quackmode = 'beg',
        quackincrement = False
        )

    flagdata(
        vis = basename+'-selected.ms',
        mode = 'manualflag',
        selectdata = True,
        antenna = 'VA08, VA11'
        )

    flagdata(
        vis = basename+'-selected.ms',
        mode = 'manualflag',
        selectdata = True,
        timerange='2:48:25~mytime_here' # problem 3
        )
    
    flagdata(
        vis = basename+'-selected.ms',
        mode = 'manualflag',
        selectdata = True,
        timerange='02:50:45.0'
        )
    
    flagdata(
        vis = basename+'-selected.ms',
        mode = 'manualflag',
        selectdata = True,
        timerange='03:10:45.0'
        )
    
    flagdata(
        vis = basename+'-selected.ms',
        mode = 'manualflag',
        selectdata = True,
        antenna = '0015' 
    )

mystep = 3
if(mystep in mysteps or doall):
    print "Step ", mystep, steptitle[mystep]

    clearcal(
        vis = basename+'-selected.ms'
        )
        
    # 1331+305 a.k.a. 3C286 is a flux calibrator known to CASA
    setjy(
        vis = basename+'-selected.ms',
        field = '2',
        spw = '0,1',   # all spws
        modimage = os.environ["CASAPATH"].split()[0]+'/data/nrao/VLA/CalModels/3C286_K.im',
        # cal. source is resolved and there is an image for it in CASA
        fluxdensity = myoption_here, # problem 4
        standard = 'Perley-Taylor 99'
        )

    # determine antenna-based complex gain for the flux and phase calibrator
    #   type G, polarization-dependent

    gaincal(
        vis = basename+'-selected.ms',
        caltable = basename+'-field1+2-gaincal.G',
        field = '1,2', # both flux and phase calibrator
        spw = '', # all spws
        calmode = 'ap', # amplitude and phase solution (default)
        gaintype = 'G',
        solint = 'inf', 
        combine = '', # calibrate the two spws separately
        refant = 'VA25', # name, not index
        minsnr = 1.0, # minimum signal to noise
        gaincurve = True, # internal VLA gain curve correction
        opacity = 0.0
        )
    
    plotcal(
        caltable = basename+'-field1+2-gaincal.G',
        xaxis = 'time',
        yaxis = 'amp',
        iteration = 'antenna'
        )
    user_check = raw_input('Quit plotcal, then press return ...')
    
    plotcal(
        caltable = basename+'-field1+2-gaincal.G',
        xaxis = 'time',
        yaxis = 'phase',
        iteration = 'antenna'
        )
    user_check = raw_input('Quit plotcal, then press return ...')
    
    # transfer the flux density scale to the phase calibrator

    os.system('rm -rf '+basename+'-field1+2-gaincal.Gflx')
    fluxscale(
        vis = basename+'-selected.ms',
        caltable = basename+'-field1+2-gaincal.G',
        fluxtable = basename+'-field1+2-gaincal.Gflx',
        reference = mystring_here, # the flux cal - problem 5
        transfer = ['12191+48299'] # the phase cal
        )

    # apply it to all the data
    
    applycal(
        vis = basename+'-selected.ms',
        gaintable = basename+'-field1+2-gaincal.Gflx',
        field = '',
        spw = '',
        gaincurve = myoption_here, # problem 6
        opacity = 0.0,
        spwmap = []
        )

# clean the phase calib combining both sidebands
mystep = 4
if(mystep in mysteps or doall):
    print "Step ", mystep, steptitle[mystep]

    os.system('rm -rf '+basename+'-phasecalib*')

    clean(
        vis = basename+'-selected.ms',
        imagename= basename+'-phasecalib',
        field=mystring_here,  # problem 7
        spw='',
        cell=myvalue_here, # problem 7 (hint should be ca. 3 times smaller than the resolution = lambda/3000m)
        imsize=[128,128],
        stokes='I',
        mode='channel',nchan=1,start=0,width=2,
        outframe='',
        psfmode='clark',
        imagermode='csclean',
        scaletype='SAULT',
        niter=mynumber_here, # problem 7
        threshold='1.0mJy',
        restfreq='22235.0MHz', # the rest freq of the H2O maser line
        phasecenter='',
        weighting='natural',
        interactive=F,
        minpb=0.1,
        pbcor=F)

mystep = 5
if(mystep in mysteps or doall):
    print "Step ", mystep, steptitle[mystep]    

    # problem 8: use imhead to list image properties
    
# clean the target combining both SPWs (sidebands)
mystep = 6
if(mystep in mysteps or doall):
    print "Step ", mystep, steptitle[mystep]

    os.system('rm -rf '+basename+'-target*')

    clean(
        vis = basename+'-selected.ms',
        imagename=basename+'-target',
        field='0',
        spw='',
        cell=myvalue_here, # problem 9
        imsize=[128,128],
        stokes='I',
        mode='channel', nchan=1, start=0, width=2, 
        outframe='',
        psfmode='clark',
        imagermode='csclean',
        scaletype='SAULT',
        niter=1000,
        gain = 0.1,
        restoringbeam = mybeam_here, # problem 9
        threshold='0.5mJy',
        restfreq='22235.0MHz', # the rest freq of the H2O maser line
        phasecenter='',
        weighting='natural',
        interactive=F,
        minpb=0.1,
        pbcor=T
        )

mystep = 7
if(mystep in mysteps or doall):
    print "Step ", mystep, steptitle[mystep]    

    # problem 10: add call to imfit here