##########################################################################
# 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('reduce-m100-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('reduce-m100-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  0: Write a flagdata2 command to flag the channels 
#          239, 447/448, 720/721 and 2847/2848 in all SPWs
# Step  4: Flag the CO line of Titan and noisy data at uvdist < 40 klambda
# Step  5: Determine the solint parameter in gaincal
# Step  6: Determine the field, solint, and calmode parameters in gaincal
# Step  7: Determine the gaintable parameter in gaincal
# Step  8: Determine the right values for reference and transfer in fluxscale
# Step  9: Complete the last applycal command
# Step 13: Determine the missing parameters in split
# Step 15: Determine the solint parameter in gaincal
# Step 17: Determine the mode, imagermode, and mask parameters in clean
# Step 18: Determine the fitspw parameter in uvcontsub
# Step 19: The mask test-M100line-orig.mask is missing. Generate it using
#          interactive clean.
# Steps 19+20: Determine the field, outframe, and weighting parameters in clean
# Step 21: Determine the axis and includepix parameters in immoments
#
##########################################################################

# EU ALMA Regional Centre CASA Tutorial for ALMA Cycle 0, January 2012
# Sample ALMA data analysis script
# "M100" 

step_title = { 0: 'Flagging \n'\
               ' (Certain channels are bad, certain antennas have problems during a few short periods)',
	       1: 'Rebin to a reduced resolution of approx. 10 km/s \n'\
               ' (in order to save processing time in the following steps)',
	       2: 'Fast phase-only gaincal for bandpass \n'\
               ' (This is essentially a one-step self-cal on the bandpass calibrator to take out fast phase fluctuations)',
               3: 'Bandpass \n'\
               ' (Frequency-dependence calibration: applying the phase cal table from the previous step on-the-fly, \n'\
               ' we generate the bandpass cal table)',
               4: 'Setjy \n'\
               ' (Titan is not a point source. After flagging a strong CO line in channels 200~479 of Titan, we set \n'\
               ' the correct model for Titan in all SPWs)',
               5: 'Fast phase-only gaincal \n'\
               ' (Similar to step 2, we do a phase-only selfcal on the bandpass, the primary phase, and the flux \n'\
               ' calibrator to improve S/N)',
               6: 'Slow phase-only gaincal \n'\
               ' (The phase-only cal to be used for the secondary phase calibrator and the target. Long solint in order \n'\
               ' to improve stability of the interpol.solution for the target)',
               7: 'Slow amp and phase gaincal \n'\
               ' (Time-dependence calibration: using the bandpass solution from step 3 and the fast phase-only solution \n'\
               '  from step 5 on-the-fly)',
               8: 'Fluxscale \n'\
               ' (using the amp and phase solution from step 7 on-the-fly, the flux scale is transferred from Titan to \n'\
               ' the phase and the bandpass calibrator)',
               9: 'Applycal \n'\
               ' (apply the bandpass solution from step 3, the phase solutions from steps 5 or 6 respectively, and the \n'\
               ' flux solution from step 8)',
               10: 'Test image of the secondary phase cal',
               11: 'Test image of the primary phase cal',
               12: 'Test image of Titan',
               13: 'Split out corrected data and time average \n'\
               ' (Now that the time dependence is calibrated, we can reduce the time resolution to 60s and save time\n'\
               ' in the following steps)',
               14: 'Concatenate data \n'\
               ' (To prepare imaging, the two datasets X220 and X54 are merged)',
               15: 'Adjust fluxscale \n'\
               ' (Since we calibrated the flux scale independently for X220 and X54, we need to align the two flux \n'\
               ' calibrations assuming that the flux of Titan did not change within the short time between X220 and X54)',
               16: 'Split out the corrected M100 data \n'\
               ' (in order to further shring the input MS for imaging)',
               17: 'Continuum image of M100 \n'\
               ' (excluding the channels of the CO(1-0) line determined with plotms)',
               18: 'Determine and subtract continuum \n'\
               ' (only for spw 0 where the CO line is)',
               19: 'Test image of central field \n'\
               ' (averaging over the CO line channels)',
               20: 'Clean CO(1-0) line cube mosaic \n'\
               ' (for the CO line channels)',
               21: 'Make moment maps \n'\
               ' (for the integrated sum and the gradient - the velocity field - of the cube;\n'\
               ' also demonstrate how to make overlay plots of the maps using imview)'
               }

# global defs
basename=['X54','X220']
makeplots=True
therefant = 'DV01'
    
############################################################################

# 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.'
if(thesteps==[]):
    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 
            
# flagging
mystep = 0
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:

        flagmanager(vis=name+'-line.ms', mode='restore', versionname='apriori')

        # Edge channels
        flagdata2(vis=name+'-line.ms', selectdata=T, 
                  field='',  manualflag=T, 
                  mf_spw='0~3:0~10;3800~3839',
                  flagbackup = F)

        # Channels 239, 447/448, 720/721 and 2847/2848 are off in all SPWs
        # PROBLEM: write a flagdata2 command to flag these channels

    # some ints are off
    flagdata2(
	    vis='X220-line.ms',
	    selectdata=F,
	    manualflag=T,
	    mf_timerange='19:52:55~19:53:04',
	    flagbackup = F)
    flagdata2(
	    vis='X54-line.ms',
	    mf_antenna='PM01',
	    mf_timerange='19:03:35~19:03:42',
	    selectdata=F,
	    manualflag=T,
	    flagbackup = F)

    flagdata2(
	    vis='X54-line.ms',
	    mf_antenna='DV04',
	    mf_timerange='19:38:45~19:38:55',
	    selectdata=F,
	    manualflag=T,
	    flagbackup = F)
    
    timing()

# Bin it up to lower spectral resolution to about 10 km/s
mystep = 1
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf '+name+'-line-vs.ms')
        split(vis=name+'-line.ms',
              outputvis=name+'-line-vs.ms',
              datacolumn='data',
              width='8')

        flagmanager(vis=name+'-line-vs.ms', mode='save', versionname='apriori')

    timing()


# Fast phase-only gaincal for bandpass
mystep = 2
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf cal-'+name+'-BPint.Gp')
        gaincal(vis=name+'-line-vs.ms', 
                caltable='cal-'+name+'-BPint.Gp', 
                spw='*:190~310',
                field='*Bandpass*',
                selectdata=F, solint='int', refant=therefant, calmode='p')

        if(makeplots):
            plotcal(caltable='cal-'+name+'-BPint.Gp', 
                    xaxis = 'time', yaxis = 'phase',
                    poln='X', plotsymbol='o', plotrange = [0,0,-180,180], 
                    iteration = 'spw',
                    figfile='cal-'+name+'-phase_vs_time_XX.BPint.Gp.png', subplot = 221)
		
            plotcal(caltable='cal-'+name+'-BPint.Gp', 
                    xaxis = 'time', yaxis = 'phase',
                    poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], 
                    iteration = 'spw',
                    figfile='cal-'+name+'-phase_vs_time_YY.BPint.Gp.png', subplot = 221)

    timing()

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

    for name in basename:
        os.system('rm -rf cal-'+name+'.B1')
        
        bandpass(vis=name+'-line-vs.ms', 
                 caltable='cal-'+name+'.B1', 
                 field='*Bandpass*',
                 bandtype='B', fillgaps=10, solnorm = T, combine='',
                 selectdata=F,
                 solint='inf',
                 refant=therefant,
                 gaintable='cal-'+name+'-BPint.Gp')
	
        if(makeplots):
            for spw in ['0','1','2','3']:
                plotcal(caltable = 'cal-'+name+'.B1', 
                        xaxis='freq', yaxis='phase', spw=spw, antenna='',
                        iteration='antenna',
                        subplot=431, overplot=False, plotrange = [0,0,-70,70],
                        plotsymbol='.', timerange='',
                        figfile='cal-'+name+'-phase.spw'+spw+'.B1.png')

                plotcal(caltable = 'cal-'+name+'.B1', 
                        xaxis='freq', yaxis='amp', spw=spw,
                        iteration='antenna',
                        subplot=431, overplot=False,
                        plotsymbol='.', timerange='',
                        figfile='cal-'+name+'-amplitude.spw'+spw+'.B1.png')
    timing()
	

# Setjy
# Strong line for Titan is obvious
# Noisy for uvdistances less than 40 klambda
mystep = 4
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        flagmanager(vis=name+'-line-vs.ms', mode='restore', versionname='apriori')

	flagdata2(vis=name+'-line-vs.ms',
                  selectdata=T, 
		  field=PROBLEM,  manualflag=PROBLEM, 
		  mf_uvrange=[PROBLEM,PROBLEM], 
		  mf_spw=[PROBLEM,PROBLEM],
		  flagbackup = F)
 
	setjy(vis=name+'-line-vs.ms',
              field='Titan',
              standard='Butler-JPL-Horizons 2010', 
	      spw='0,1,2,3')

    timing()

# Fast phase-only gaincal 
mystep = 5
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf cal-'+name+'-int.Gp')
	gaincal(vis=name+'-line-vs.ms', 
		caltable='cal-'+name+'-int.Gp', 
		spw='*:25~455', 
		field='*Phase*,*Band*,Titan',
		gaintable='cal-'+name+'.B1',
		selectdata=F, solint='PROBLEM', 
		refant=therefant, calmode='p')

	if(makeplots):
            plotcal(caltable='cal-'+name+'-int.Gp', 
                    xaxis = 'time', yaxis = 'phase',
                    poln='X', plotsymbol='o', plotrange = [0,0,-180,180], 
                    iteration = 'spw',
                    figfile='cal-'+name+'-phase_vs_time_XX.int.Gp.png', subplot = 221)

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

    timing()


# Slow phase-only gaincal 
mystep = 6
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf cal-'+name+'-scan.Gp')
	gaincal(vis=name+'-line-vs.ms', 
		caltable='cal-'+name+'-scan.Gp', 
		spw='*:25~455', 
		field='PROBLEM',
		gaintable='cal-'+name+'.B1',
		selectdata=F, solint='PROBLEM', 
		refant=therefant, calmode='PROBLEM')

	if(makeplots):
		plotcal(caltable='cal-'+name+'-scan.Gp', 
			xaxis = 'time', yaxis = 'phase',
			poln='X', plotsymbol='o', plotrange = [0,0,-180,180], 
			iteration = 'spw',
			figfile='cal-'+name+'-phase_vs_time_XX.scan.Gp.png', subplot = 221)

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

    timing()


# Slow amplitude and phase gaincal 
mystep = 7
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf cal-'+name+'-scan.Gap')
	gaincal(vis=name+'-line-vs.ms', 
		caltable='cal-'+name+'-scan.Gap', 
		spw='*:25~455', 
		field='*Phase*,*Band*,Titan',
		gaintable=['PROBLEM'],
		selectdata=F, solint='inf', 
		refant=therefant, calmode='ap')

	if(makeplots):
		plotcal(caltable='cal-'+name+'-scan.Gap', 
			xaxis = 'time', yaxis = 'amp',
			poln='X', plotsymbol='o', 
			iteration = 'spw',
			figfile='cal-'+name+'-amp_vs_time_XX.scan.Gap.png', subplot = 221)

		plotcal(caltable='cal-'+name+'-scan.Gap', 
			xaxis = 'time', yaxis = 'amp',
			poln='Y', plotsymbol='o', 
			iteration = 'spw',
			figfile='cal-'+name+'-amp_vs_time_YY.scan.Gap.png', subplot = 221)

    timing()

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

    for name in basename:
	fluxscale(vis=name+'-line-vs.ms',caltable='cal-'+name+'-scan.Gap',
		  fluxtable='cal-'+name+'.flux',reference='PROBLEM',transfer='PROBLEM')

    timing()


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

    for name in basename:
	# to the bandpass cal
	applycal(vis=name+'-line-vs.ms',field='*Band*',
	        gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp','cal-'+name+'.flux'],
	        interp=['nearest','nearest','nearest'],
	        gainfield=['*Band*','*Band*','*Band*'],
		calwt=F,
		flagbackup=T)

	# to the secondary phase cal
	applycal(vis=name+'-line-vs.ms',field='3c273 - Phase',
	        gaintable=['cal-'+name+'.B1','cal-'+name+'-scan.Gp','cal-'+name+'.flux'],
       	 	interp=['nearest','linear','linear'],
       	 	gainfield=['*Band*','1224*','1224*'],
		calwt=F,
		flagbackup=T)

	# to the primary phase cal
	applycal(vis=name+'-line-vs.ms',field='1224*',
        	gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp','cal-'+name+'.flux'],
        	interp=['nearest','nearest','nearest'],
        	gainfield=['*Band*','1224*','1224*'],
		calwt=F,
		flagbackup=T)

	# to Titan
        applycal(vis=name+'-line-vs.ms',field='Titan',
                 gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp','cal-'+name+'.flux'],
                 interp=['nearest','nearest','nearest'],
                 gainfield=['*Band*','Titan','Titan'],
                 calwt=F,
                 flagbackup=T)

	# to M100
        applycal(PROBLEM)
        

    timing()


# Test image of the secondary phase cal
mystep = 10
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf test-'+name+'-sec_phasecal*')
	clean(vis=name+'-line-vs.ms',
	      imagename='test-'+name+'-sec_phasecal',
	      field='3c*Ph*',spw='0~3',
	      nterms=2,
	      mode='mfs',niter=100,
	      interactive=False,
	      mask=[96, 96, 104, 104],
              imsize=200,cell='0.5arcsec')
    if makeplots:
        for name in basename:
            imview(raster={'file': 'test-'+name+'-sec_phasecal.image.tt0', 'colorwedge':T,
                           'range':[-0.02, 8.0], 'scaling':-1.5, 'colormap':'Rainbow 2'},
                   out='test-'+name+'-sec_phasecal.png', zoom=1)

    timing()


# Test image of the primary phase cal
mystep = 11
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf test-'+name+'-prim_phasecal*')
	clean(vis=name+'-line-vs.ms',
	      imagename='test-'+name+'-prim_phasecal',
	      field='1224*',spw='0~3',
	      nterms=2,
	      mode='mfs',niter=100,
	      interactive=False,
	      mask=[96, 96, 104, 104],
              imsize=200,cell='0.5arcsec')

    if makeplots:
        for name in basename:
            imview(raster={'file': 'test-'+name+'-prim_phasecal.image.tt0', 'colorwedge':T,
                           'range':[-0.005, 0.9], 'scaling':-2.5, 'colormap':'Rainbow 2'},
                   out='test-'+name+'-prim_phasecal.png', zoom=1)
            calstat=imstat(imagename='test-'+name+'-prim_phasecal.image.tt0', region='', box='30,30,170,80')
            rms=(calstat['rms'][0])
            print '>> rms in phase calibrator image: '+str(rms)
            calstat=imstat(imagename='test-'+name+'-prim_phasecal.image.tt0', region='')
            peak=(calstat['max'][0])
            print '>> Peak in phase calibrator image: '+str(peak)
            print '>> Dynamic range in phase calibrator image: '+str(peak/rms)

    timing()

# Test image of Titan
mystep = 12
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf test-'+name+'-Titan*')
	clean(vis=name+'-line-vs.ms',
	      imagename='test-'+name+'-Titan',
	      field='Titan',spw='0~3',
	      mode='mfs',niter=100,
	      interactive=False,
	      mask=[96, 96, 104, 104],imsize=200,cell='0.5arcsec')

    timing()

# Split out corrected data and time average at the same time into 1 minute time bins
mystep = 13
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf '+name+'-corrected.ms*')
        split(vis=name+'-line-vs.ms', outputvis=name+'-corrected.ms'
              ) # PROBLEM

    timing()
    

# Concat
mystep = 14
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf all.ms*')
    concat(vis=['X54-corrected.ms', 'X220-corrected.ms'],
           concatvis='all.ms', copypointing=False)

    timing()

# Adjust fluxscale
mystep = 15
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# INFO	fluxscale::::	 Flux density for 1224+213 Phase in SpW=0 is: 1.07959 +/- 0.0105176 (SNR = 102.646, N= 24)
# INFO	fluxscale::::	 Flux density for 1224+213 Phase in SpW=1 is: 1.09821 +/- 0.0102055 (SNR = 107.61, N= 24)
# INFO	fluxscale::::	 Flux density for 1224+213 Phase in SpW=2 is: 1.11718 +/- 0.00845822 (SNR = 132.082, N= 24)
# INFO	fluxscale::::	 Flux density for 1224+213 Phase in SpW=3 is: 1.12509 +/- 0.00724378 (SNR = 155.318, N= 24)

# INFO    listobs::ms::summary+     SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)  TotBW(kHz)  Corrs
# INFO    listobs::ms::summary+     0         480 TOPO  113728.128  3906.25       1875000     XX  YY
# INFO    listobs::ms::summary+     1         480 TOPO  111853.128  3906.25       1875000     XX  YY
# INFO    listobs::ms::summary+     2         480 TOPO  103661.722  -3906.25      1875000     XX  YY
# INFO    listobs::ms::summary+     3         480 TOPO  101849.222  -3906.25      1875000     XX  YY

    # calculate spectral index of primary phase cal 
    # (alternatively run setjy separately on each SPW)
    freq0 = 101849.222 - 480./2. * 3906.25E-3
    freq1 = 113728.128 + 480./2. * 3906.25E-3
    flux0 = 1.12509
    flux1 = 1.07959
    myspix = log(flux0/flux1)/log(freq0/freq1)
    print "Spectral index of primary phase cal: ", myspix

    setjy(vis='all.ms',
          field='1224*',
          fluxdensity=flux0,
          reffreq=str(freq0)+'MHz',
          spix=myspix,
          spw='0,1,2,3')

    os.system('rm -rf cal-fluxadjust.Ga')
    gaincal(vis='all.ms', 
            caltable='cal-fluxadjust.Ga', 
            spw='*:25~455', 
            field='1224*',
            gaintable=[],
            selectdata=F,
            solint='PROBLEM', 
            refant=therefant, calmode='a')

    if(makeplots):
        plotcal(caltable='cal-fluxadjust.Ga', 
                xaxis = 'time', yaxis = 'amp',
                poln='X', plotsymbol='o', 
                iteration = 'spw',
                figfile='cal-all-amp_vs_time_XXfluxadjust.Ga.png', subplot = 221)

    applycal(vis='all.ms',field='M100',
             gaintable=['cal-fluxadjust.Ga'],
             interp=['linear'],
             gainfield=['1224*'],
             calwt=F,
             flagbackup=T)

    timing()

# split out the corrected M100 data 
mystep = 16
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf M100all_lores.ms*')
    split(vis='all.ms', outputvis='M100all_lores.ms',
          field = 'M100', 
          datacolumn='corrected')

    timing()


# Continuum image
mystep = 17
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf M100cont.*')
    clean(vis = 'M100all_lores.ms',
          imagename = 'M100cont',
          field='2~47',
          spw='0:10~210;256~440,1~3:10~460', # exclude CO(1-0) line emission 
          mode = 'PROBLEM',
          niter = 1000,
          mask=[0,0,0,0],# PROBLEM
          imagermode = 'PROBLEM', 
          interactive = F,
          imsize = 200,
          cell = '0.5arcsec',
          phasecenter='J2000 12h22m54.9 +15d49m15')

# Continuum peak is 0.5 mJy. Too weak for self-cal...

    timing()


# uvcontsub2
mystep = 18
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]
    
    os.system('rm -rf M100all_lores.ms.c*')
    uvcontsub(vis='M100all_lores.ms',field='',fitspw='PROBLEM',
              combine='',solint='inf',fitorder=1,spw='0',want_cont=False)

    timing()


# Test image of central field
mystep = 19
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf test-M100line.*')
    clean(
	vis='M100all_lores.ms.contsub',
        imagename='test-M100line',
	field='PROBLEM',
	spw='0:231~248',
	mode='mfs',
	niter=500,gain=0.1,threshold='0.0mJy',
	imagermode='csclean',
	mask='test-M100line-orig.mask', # PROBLEM: mask missing
	interactive=False, # switch to interactive to determine it
	outframe='PROBLEM',veltype='radio',
	imsize=200,cell='0.5arcsec',
	phasecenter='',
	stokes='I',
	weighting='PROBLEM',robust=0.5,
	calready=False,
	npercycle=100,cyclefactor=1.5,cyclespeedup=-1)

    timing()


# Clean line cube mosaic
mystep = 20
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf M100line.*')
    clean(vis='M100all_lores.ms.contsub',imagename='M100line',
          field='PROBLEM',
          spw='0:220~259',
          mode='channel',
          niter=1000,gain=0.1,threshold='0.0mJy',psfmode='clark',
          imagermode='mosaic',ftmachine='mosaic',mosweight=False,
          scaletype='SAULT',
          interactive=False,
          mask='M100line-orig.mask',
          nchan=40,start=220,
          width=1,
          outframe='PROBLEM',veltype='radio',
          imsize=600,cell='0.5arcsec',
          phasecenter='J2000 12h22m54.9 +15d49m10',
          restfreq='115.271201800GHz',stokes='I',
          weighting='PROBLEM',robust=0.5,
          pbcor=False,minpb=0.2,
          calready=False,
          npercycle=100,cyclefactor=1.5,cyclespeedup=-1)

    timing()

# Moments
mystep = 21
if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf M100-CO.mom?')
    immoments(imagename='M100line.image',
              moments=[0],
              axis='PROBLEM',
              region='',box='100,110,515,500',
              chans='7~35',
              mask='',
              outfile='M100-CO.mom0',
              includepix=[PROBLEM])

    immoments(imagename='M100line.image',
              moments=[1],
              axis='PROBLEM',
              region='',box='100,110,515,500',
              chans='7~35',
              mask='',
              outfile='M100-CO.mom1',
              includepix=[PROBLEM])

    if makeplots:
        os.system('rm -rf M100-CO_velfield.png')
        imview(contour={'file': 'M100-CO.mom0','levels': 
                        [1,2,5,10,20,40,80,160],'base':0,'unit':1}, 
               raster={'file': 'M100-CO.mom1','range': [1440,1700],
                       'colorwedge':T, 'colormap': 'Rainbow 2'}, out='M100-CO_velfield.png')
        os.system('rm -rf M100-CO_map.png')
        imview(contour={'file': 'M100-CO.mom1','levels': 
                        [1430,1460,1490,1520,1550,1580,1610,1640,1670,1700],'base':0,'unit':1}, 
               raster={'file': 'M100-CO.mom0', 'colorwedge':T,
                       'colormap': 'Rainbow 2','scaling':-1.8,'range': [0.5,20]}, 
               out='M100-CO_map.png')

        os.system('rm -rf M100-CO_contmap.png')
        imview(contour={'file': 'M100cont.image','levels': 
                        [0.00025,0.0004],'base':0,'unit':1}, 
               zoom=3,
               raster={'file': 'M100-CO.mom0', 'colorwedge':T,
                       'colormap': 'Rainbow 2','scaling':0,'range': [0.8,40]}, 
               out='M100-CO_contmap.png')

    timing()

print 'Done.'