##########################################################################
# This casapy script was prepared for the EU ARC Cycle 0 PI CASA tutorial.
# It was first made to work correctly.
# Then, a number of problems were introduced.
# It is your task to fix them.
# You will find a list  of the problems further below.
# Look for the string 'PROBLEM' to find all problems. 
#
# The script uses some infrastructure to make it easier for you
# to run and repeat the script step by step:
# Further below you find the "step_title" Python dictionary which
# shows you what is happening in each analysis step.
# You can start this script by starting casapy and entering
#     execfile('selfcal-ngc3256-with-problems.py')
# By default it will run all steps.
# (Of course it will only run correctly when you have fixed all problems.)
# If you would like to run individual steps, you can set the
# variable mysteps before starting the script, e.g.:
#
#     mysteps = [4,5,6]
#     execfile('selfcal-ngc3256-with-problems.py')
#
# This example would only execute steps 4, 5, and 6.
# Setting mysteps to an empty list [] will execute all steps.
#
#
# List of problems in the individual analysis steps of this script:
# Step 1: Determine the spw, solint, refant, and calmode
# Step 2: Determine the interp and gaintable parameters
# Step 3: Determine the spw, mode, threshold, and mask parameters
# Step 4: Determine the fitspw, solint, fitorder, and combine parameters
# Step 6: Determine the outframe and resfreq parameters
# Step 7: Determine the chans and box parameters
# Step 8: Determine the chans and box parameters
#
##########################################################################


# EU ALMA Regional Centre CASA Tutorial for ALMA Cycle 0
# Sample ALMA data PI analysis script
# "NGC 3256 self-calibration"
#
# Get the input MS from http://almascience.eso.org/almadata/sciver/NGC3256/NGC3256_Band3_CalibratedData.tgz
#

step_title = { 0: 'Continuum image of the galaxy',
	       1: 'Phase-only gaincal for first iteration of selfcal',
               2: 'Applycal of phase solutions',
               3: 'Clean with new cal - get selfcalibrated cont image',
               4: 'Continuum subtraction',
               5: 'Determine restfreq for CO line in spw 0',
               6: 'Clean individual channels - get selfcalibrated line-cube for spw 0',
               7: 'Make moment maps',
               8: 'Make velocity dispersion map'
               }

# global defs
makeplots=True
    
############################################################################

# Some infrastructure to make repeating individual parts
#   of this workflow more convenient.

totaltime = 0
inittime = time.time()
ttime = inittime
steptime = []
thesteps = []

def timing():
    global totaltime
    global inittime
    global ttime
    global steptime
    global step_title
    global mystep
    global thesteps
    thetime = time.time()
    dtime = thetime - ttime
    steptime.append(dtime)
    totaltime += dtime
    ttime = thetime
    casalog.origin('TIMING')
    casalog.post( 'Step '+str(mystep)+': '+step_title[mystep], 'WARN')
    casalog.post( 'Time now: '+str(ttime), 'WARN')
    casalog.post( 'Time used this step: '+str(dtime), 'WARN')
    casalog.post( 'Total time used so far: ' + str(totaltime), 'WARN')
    casalog.post( 'Step  Time used (s)  Fraction of total time (percent) [description]', 'WARN')
    for i in range(0, len(steptime)):
        casalog.post( '  '+str(thesteps[i])+'   '+str(steptime[i])+'  '+str(steptime[i]/totaltime*100.)
                      +' ['+step_title[thesteps[i]]+']', 'WARN')
        
def asyncwait(handles=[]):
    done = (handles==[])
    while not done:
        for handle in handles:
            done = (tm.retrieve(handle)['status'] == 'done')
            if not done:
                break
        time.sleep(10)
    return


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

# The Python variable 'mysteps' will control which steps
# are executed when you start the script using
#   execfile('ngc4258_tutorial.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.

##########################################################################

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

    if makeplots:

        plotms(vis='ngc3256_line_target.ms',spw='0',xaxis='channel',yaxis='amp',
               avgtime='1e8',avgscan=T, plotfile='plot-amp-vs-chan-co-line.png',
               overwrite=True)
        plotms(vis='ngc3256_line_target.ms',spw='1',xaxis='channel',yaxis='amp',
               avgtime='1e8',avgscan=T, plotfile='plot-amp-vs-chan-cn-lines.png',
               overwrite=True)
        plotms(vis='ngc3256_line_target.ms',spw='2',xaxis='channel',yaxis='amp',
               avgtime='1e8',avgscan=T, plotfile='plot-amp-vs-chan-linefree2.png',
               overwrite=True)
        plotms(vis='ngc3256_line_target.ms',spw='3',xaxis='channel',yaxis='amp',
               avgtime='1e8',avgscan=T, plotfile='plot-amp-vs-chan-linefree3.png',
               overwrite=True)

        # If you would like the script to stop at this point, you could insert the following line
        # user_check=raw_input('press enter to continue script\n')

    os.system('rm -rf result-ngc3256_cont*')
    clean( vis='ngc3256_line_target.ms', imagename='result-ngc3256_cont',
           spw='0:20~53;71~120,1:70~120,2:20~120,3:20~120', psfmode='hogbom',
           mode='mfs', niter=500, threshold='0.3mJy',imsize=100, cell='1arcsec',
           weighting='briggs', robust=0.0, interactive=T)

#   clean: PROBLEM (optional) save mask for future use

    calstat=imstat(imagename='result-ngc3256_cont.image', region='', box='10,10,90,35')
    rms=(calstat['rms'][0])
    print '>> rms in continuum image: '+str(rms)
    calstat=imstat(imagename='result-ngc3256_cont.image', region='')
    peak=(calstat['max'][0])
    print '>> Peak in continuum image: '+str(peak)
    print '>> Dynamic range in continuum image: '+str(peak/rms)

    if makeplots:

        imview(raster={'file': 'result-ngc3256_cont.image', 'colorwedge':T,
                       'range':[-0.001, 0.009], 'scaling':0, 'colormap':'Rainbow 2'},
               out='result-ngc3256_cont.png', zoom=2)

    timing()


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

    os.system('rm -rf cal-ngc3256_cont_30m.Gp')

    gaincal(vis='ngc3256_line_target.ms', field='NGC*',
            caltable='cal-ngc3256_cont_30m.Gp',
            spw='PROBLEM', solint='PROBLEM', combine='scan',
            refant='PROBLEM', calmode='PROBLEM',
            minblperant=3)

    if makeplots:

        plotcal(caltable = 'cal-ngc3256_cont_30m.Gp', xaxis = 'time', yaxis = 'phase',
                poln='X', plotsymbol='o', plotrange = [0,0,-180,180],
                iteration = 'spw', figfile='cal-phase_vs_time_XX_30_Gp.png',
                subplot = 221)

        plotcal(caltable = 'cal-ngc3256_cont_30m.Gp', xaxis = 'time', yaxis = 'phase',
                poln='Y', plotsymbol='o', plotrange = [0,0,-180,180],
                iteration = 'spw', figfile='cal-phase_vs_time_YY_30_Gp.png',
                subplot = 221)

    timing()

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

    applycal(vis='ngc3256_line_target.ms', interp='PROBLEM',
             gaintable='PROBLEM')

    timing()

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

    os.system('rm -rf result-ngc3256_cont_sc1*')
    clean( vis='ngc3256_line_target.ms', imagename='result-ngc3256_cont_sc1',
           spw='PROBLEM', psfmode='hogbom', mode='PROBLEM', niter=500,
           threshold='PROBLEM', mask='PROBLEM', imsize=100, cell='1arcsec',
           weighting='briggs', robust=0.0, interactive=T)

    calstat=imstat(imagename='result-ngc3256_cont_sc1.image', region='', 
                   box='10,10,90,35')
    rms=(calstat['rms'][0])
    print '>> rms in continuum image: '+str(rms)
    calstat=imstat(imagename='result-ngc3256_cont_sc1.image', region='')
    peak=(calstat['max'][0])
    print '>> Peak in continuum image: '+str(peak)
    print '>> Dynamic range in continuum image: '+str(peak/rms)

    if makeplots:
        imview(raster={'file': 'result-ngc3256_cont_sc1.image', 'colorwedge':T,
                       'range':[-0.001, 0.009], 'scaling':0, 'colormap':'Rainbow 2'},
               out='result-ngc3256_cont_sc1.png', zoom=2)

    timing()


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

    os.system('rm -rf ngc3256_line_target.ms.contsub')
    uvcontsub(vis = 'ngc3256_line_target.ms',
              fitspw='PROBLEM', solint ='PROBLEM', 
              fitorder = 'PROBLEM',
              combine='PROBLEM')
    
    timing()

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

    os.system('rm -rf myresults.tbl')
    slsearch(outfile='myresults.tbl', freqrange = [84,116], species=['COv=0'])
    sl.open('myresults.tbl')
    sl.list()

    tb.open('myresults.tbl') 
    restfreq=tb.getcol('FREQUENCY')[0]
    tb.close()
    print 'COv=0 line rest frequency = ', restfreq, ' GHz'

    timing()

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

    os.system('rm -rf result-ngc3256_line_CO.*')
    clean(vis='ngc3256_line_target.ms.contsub', imagename='result-ngc3256_line_CO',
          spw='0:38~87', mode='channel', start='', nchan=50, width='', 
          psfmode='hogbom', outframe='PROBLEM', restfreq='PROBLEM', 
          mask=[53,50,87,83], niter=2000, interactive=F, imsize=128, cell='1arcsec', 
          weighting='briggs', robust=0.0, threshold='5mJy')

    timing()


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

    os.system('rm -rf result-ngc3256_CO1-0.mom0*')
    immoments(imagename='result-ngc3256_line_CO.image', moments=[0],
              chans='PROBLEM', box='PROBLEM', axis='spectral',
              includepix=[0.02, 10000], outfile='result-ngc3256_CO1-0.mom0')

    os.system('rm -rf result-ngc3256_CO1-0.mom1*')
    immoments(imagename='result-ngc3256_line_CO.image', moments=[1],
              chans='PROBLEM', box='PROBLEM', axis='spectral',
              includepix=[0.045, 10000], outfile='result-ngc3256_CO1-0.mom1')

    if makeplots:
        imview(contour={'file': 'result-ngc3256_CO1-0.mom0','levels': 
                        [5,10,20,40,80,160],'base':0,'unit':1}, 
               raster={'file': 'result-ngc3256_CO1-0.mom1','range': [2630,2920],
                       'colorwedge':T, 'colormap': 'Rainbow 2'},
               out='result-CO_velfield.png')

        imview(contour={'file': 'result-ngc3256_CO1-0.mom1','levels': 
                        [2650,2700,2750,2800,2850,2900],'base':0,'unit':1}, 
               raster={'file': 'result-ngc3256_CO1-0.mom0', 'colorwedge':T,
                       'colormap': 'Rainbow 2','scaling':-1.0,'range': [0.8,250]}, 
               out='result-CO_map.png')

    timing()

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

    os.system('rm -rf result-ngc3256_CO1-0.mom2*')
    immoments(imagename='result-ngc3256_line_CO.image', moments=[2],
              chans='PROBLEM', box='PROBLEM', axis='spectral',
              includepix=[0.035, 10000], outfile='result-ngc3256_CO1-0.mom2')

    if makeplots:

        imview(contour={'file': 'result-ngc3256_CO1-0.mom2','levels': 
                        [20,30,40,50,60],'base':0,'unit':1}, 
               raster={'file': 'result-ngc3256_CO1-0.mom2', 'colorwedge':T,
                       'colormap': 'Greyscale 1','scaling':-1.0,'range': [0,74]}, 
               out='result-CO_dispersion.png')

    timing()