# Calibration script for NGC 3256 SciVer ALMA data in CASA 4.7.
#
# Adopted/updated for CASA 4.7 from CASA ALMA Guide at
# https://casaguides.nrao.edu/index.php/NGC3256_Band3_Calibration_for_CASA_4.3
# 
# Only final 'executive' CASA tasks are retained in the production-run script 
# - all the on-screen diagnostics and 'trial-error' steps are commented out.
#
# 2017/05/23, barta@asu.cas.cz, EU ARC, Czech node

#--- Create a list of UIDs
basename=['uid___A002_X1d54a1_X5','uid___A002_X1d54a1_X174',
          'uid___A002_X1d54a1_X2e3','uid___A002_X1d5a20_X5',
          'uid___A002_X1d5a20_X174','uid___A002_X1d5a20_X330']


#--- Print observation summary
for name in basename:
    listobs(vis=name+'.ms')

#--- Show antenna possitions
plotants(vis=basename[0]+'.ms', figfile=basename[0]+'_plotants.png')


#--- Create WVR and Tsys caltables. Remove the old ones before, for sure.
for name in basename:
    os.system('rm -rf cal-tsys_'+name+'.calnew cal-*'+name+'.Wnew')
# TSys
    gencal(vis=name+'.ms',
           caltype='tsys',
           caltable='cal-tsys_'+name+'.calnew')
# WVR
    wvrgcal(vis=name+'.ms', caltable='cal-'+name+'.Wnew',  
            toffset=-1, segsource=True,
            tie=["Titan,1037-295,NGC3256"], statsource="1037-295",
            smooth='2.88s')


#--- Apriori flagging

for name in basename:
    flagdata(vis=name+'.ms', flagbackup=False, mode = 'shadow')
    flagdata(vis=name+'.ms',mode='manual', autocorr=True)
    flagdata(vis=name+'.ms', mode='manual', flagbackup=False,
             intent='*POINTING*')
    flagdata(vis=name+'.ms', mode='manual', flagbackup=False, 
             intent='*ATMOSPHERE*')
    
# Keep the flagging done so far
    flagmanager(vis = name+'.ms', mode = 'save', versionname = 'Apriori')


#--- Inspect the data - see large phase delays (wraps) for antenna DV07.
#    These phase distlrsions are Single Band Delays (SBD) for antenna DV07.
#    This holds just for first 3 UIDs taken on April 16 2011, the next day 
#    the antenna was OK. 

#=== Commented out in production-run script

'''
plotms(vis=basename[0]+'.ms', xaxis='channel', yaxis='phase',
       spw='3', antenna='PM03&DV07', correlation='XX', avgtime='1e8', 
       showgui=True)
'''

#--- Make correction calibration table for Apr 16 data
for i in range(3): # loop over the first three ms's
    name=basename[i]
    os.system('rm -rf cal-'+name+'_del.K')
# Caltable type SBD (single-band delays)
    gencal(vis=name+'.ms', caltable='cal-'+name+'_del.Knew',
           caltype='sbd', antenna='DV07', pol='X,Y', spw='1,3,5,7',
           parameter=[1.00, 1.10, -3.0, -3.0, -3.05, -3.05, -3.05, -3.05])


#--- Apply the corrections: SBD delays & WVR phase distorsion for Apr 16 data,
#    just the WVR corrections for Apr 17 data.

for i in range(3): # loop over the first three data sets
    name=basename[i]
    applycal(vis=name+'.ms', flagbackup=False, spw='1,3,5,7',
             interp=['nearest','nearest'], 
             gaintable=['cal-'+name+'_del.Knew', 'cal-'+name+'.Wnew'])



for i in range(3,6): # loop over the last three data sets
    name=basename[i]
    applycal(vis=name+'.ms', flagbackup=False, spw='1,3,5,7',
             interp='nearest', gaintable='cal-'+name+'.Wnew')

#--- See the corrections after calibration application: Switch between "DATA"
#    and "CORRECTED" column

#=== Commented out in production-run script

'''
plotms(vis=basename[0]+'.ms', xaxis='channel', yaxis='phase',
       spw='1', antenna='', correlation='XX', avgtime='1e8',
       coloraxis='baseline', avgscan=True, selectdata=True, field='1037*', 
       showgui=True)
'''

#--- Extract the corrected data to a new MS

for name in basename:
    os.system('rm -rf '+name+'_K_WVR.ms*')
    split(vis=name+'.ms', outputvis=name+'_K_WVR.ms',
          datacolumn='corrected', spw='0~7')


#--- Inspection of Tsys data (let it run as the output is to PNG files)
for spw in ['1','3','5','7']:
    for name in basename:
        plotbandpass(caltable='cal-tsys_'+name+'.calnew', 
                     xaxis='freq', yaxis='amp',
                     spw=spw, overlay='time', # also try overlay='antenna'
                     plotrange=[0, 0, 40, 180],
                     figfile='cal-tsys_per_spw_'+spw+'_'+name+'.png', 
                     interactive=False)

#--- Scans 4,5 and 9 at antenna  DV04 in SPW7 are affected, let us flag them

flagdata(vis='uid___A002_X1d54a1_X174_K_WVR.ms', mode='manual',
         antenna='DV04', flagbackup=True, scan='4,5,9', spw='7')




#--- Apply the TSys calibration, extra for each field
for name in basename:
	for field in ['Titan','1037*','NGC*']:
		applycal(vis=name+'_K_WVR.ms', spw='1,3,5,7', 
                         flagbackup=False, field=field, gainfield=field,
                         interp='linear,spline', 
                         gaintable=['cal-tsys_'+name+'.calnew'])


#--- Extract corrected science data to separate MS - leave the not usefull
#    WVR and TSys calibration data behind
for name in basename:
    os.system('rm -rf '+name+'_line.ms*')
    split(vis=name+'_K_WVR.ms', outputvis=name+'_line.ms',
          datacolumn='corrected', spw='1,3,5,7')



#--- Concatenate the data in the one big MS

# Prepare the list of corrected partial MSs
comvis=[]
for name in basename:
	comvis.append(name+'_line.ms')
# Just for sure remove (if it exists)
os.system('rm -rf ngc3256_line.ms*')

# Make a new cocatenated MS
concat(vis=comvis, concatvis='ngc3256_line.ms')

#--- Data inspection
listobs('ngc3256_line.ms')

# Show spectrum - see the edge channels

#=== Commented out in production-run script

'''
plotms(vis='ngc3256_line.ms', xaxis='channel', yaxis='amp',
       averagedata=True, avgbaseline=True, avgtime='1e8', 
       avgscan=True, showgui=True)
'''

#--- Remove the edge channels
flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='*:0~16,*:125~127')

#=== Commented out in production-run script

# Titan has an increase of intesity becoause blending with Saturn rings
# on the next day...

'''
plotms(vis='ngc3256_line.ms', xaxis='time', yaxis='amp',
       averagedata=True, avgchannel='128', coloraxis='field',
       iteraxis='spw', showgui=True)

'''
# ...this can be seen here by source size increase.

'''
plotms(vis='ngc3256_line.ms', xaxis='uvdist', yaxis='amp', field='Titan',
       averagedata=T, avgchannel='128', avgtime='1e8', coloraxis='scan')
'''
#=== /End of diagnostic plots


#--- Flag the Saturn-ring blended data
flagdata(vis = 'ngc3256_line.ms', flagbackup=True,
	timerange='>2011/04/16/12:00:00', field='Titan')

#--- Give the Titan an ephemeris - not needed anymore as new ASDM format keeps
#    it internally 
fixplanets(vis='ngc3256_line.ms', field='Titan', fixuvw=True)



#--- Further data inspection plotms():
#   * Baselines with DV07 have very high amplitudes in spw 3, correlation YY
#   * Baselines with DV08 have very low amplitudes in spw 3, correlation YY, but
#     only for the last observation
#   * Baselines with PM03 have low amplitudes at 2011/04/17/02:15:00 for spw 0
#   * Baselines with PM03 have low amplitudes at 2011/04/16/04:15:15 for spw 2 
#     and 3
#   * For the second day, the baseline DV10-PM03 is corrupted
#
#   We flag the corresponding data

flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='3',
         correlation='YY', mode='manual',
         antenna='DV07', timerange='')

flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='3',
         correlation='YY', mode='manual',
         antenna='DV08', timerange='>2011/04/17/03:00:00')

flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='0',
         correlation='', mode='manual',
         antenna='PM03', timerange='2011/04/17/02:15:00~02:15:50')

flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='2,3', 
         correlation='', mode='manual',
         antenna='PM03', timerange='2011/04/16/04:13:50~04:18:00')

flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='',
         correlation='', mode='manual',
         antenna='PM03&DV10', timerange='>2011/04/16/15:00:00')


#--- Inspect bandpass (B) calibrator

#=== Commented out in production-run script

'''
plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='phase', selectdata=True,
       field='1037*', avgtime='1E6', avgscan=True, 
       coloraxis='baseline', iteraxis='antenna', showgui=True)

plotms(vis='ngc3256_line.ms', xaxis='time', yaxis='phase', selectdata=True,
       field='1037*', spw='0:30~90', avgchannel='128', avgscan=True, 
       coloraxis='baseline', iteraxis='antenna', showgui=True)

'''


#+++++++++++++ Inspection of primary phase (P) calibrator (1037-295) after
# application of calibration tables shows outliers fo some scans and particular
# antennas. With this apriori knowledge (found by trial-error, of course) we
# can flag the data right now to avoid repeated calibration steps.


#--- Flag corrupted data

flagdata(vis='ngc3256_line.ms', mode='manual',
	timerange='2011/04/16/04:13:35~04:13:45', flagbackup=True)

# tget(): 
# copy-paste last args of flagdata() call, we will change just some of them

tget(flagdata)


timerange='2011/04/16/05:21:13~05:21:19'
flagdata()
timerange='2011/04/16/04:16:40~04:16:49'
flagdata()
timerange='2011/04/16/04:14:00~04:17:10'; antenna='PM03'
flagdata()
timerange='2011/04/17/01:52:20~01:53:10'; antenna='DV10'
flagdata()
timerange='2011/04/17/00:35:30~01:20:20'; antenna='DV04'; spw='3'
flagdata()



#--- Short time-scale balancing of phase variations for the B calibrator

gaincal(vis='ngc3256_line.ms', caltable='cal-ngc3256.G1', 
        spw='*:40~80', field='1037*',
        selectdata=True, solint='int', refant='DV04', calmode='p')


#--- Inspect the phase caltables for B calibrator

#=== Commented out in production-run script

'''
plotcal(caltable = 'cal-ngc3256.G1', xaxis = 'time', yaxis = 'phase',
	poln='X', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw',
	figfile='cal-phase_vs_time_XX.G1.png', subplot = 221, showgui=True)


plotcal(caltable = 'cal-ngc3256.G1', xaxis = 'time', yaxis = 'phase',
	poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw',
	figfile='cal-phase_vs_time_YY.G1.png', subplot = 221, showgui=True)
'''

#--- Do bandpass calibration ('spectral balancing') separatelly for Apr 16 
#    and 17 data. Make phase corrections for B calibrator 'on the fly'
#    (instead of calling applycal() beforehand) - use parameter 'gaintable'

# Apr 16 data
bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1', 
         gaintable = 'cal-ngc3256.G1', timerange='<2011/04/16/15:00:00',
         field = '1037*', minblperant=3, minsnr=2, solint='inf', 
         combine='scan,obs',
         bandtype='B', fillgaps=1, refant = 'DV04', solnorm=True)

# Apr 17 data - append to the caltable made in previous step 
bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1', 
         gaintable = 'cal-ngc3256.G1', timerange='>2011/04/16/15:00:00',
         field = '1037*', minblperant=3, minsnr=2, solint='inf', 
         combine='scan,obs',
         bandtype='B', fillgaps=1, refant = 'DV04', solnorm=True, append=True)

#--- Plot bandpass solutions: Phase & Amplitude

#=== Commented out in production-run script
    
'''
plotbandpass(caltable = 'cal-ngc3256.B1', xaxis='freq', yaxis='phase', 
             plotrange = [0,0,-70,70],
             overlay='antenna', figfile='bandpass.B1.png')

plotbandpass(caltable = 'cal-ngc3256.B1', xaxis='freq', yaxis='amp', 
             overlay='antenna',
             figfile='bandpass.B2.png')
'''

#--- Absolute flux calibration for Titan
setjy(vis='ngc3256_line.ms', field='Titan', 
      standard='Butler-JPL-Horizons 2012', 
      spw='0,1,2,3', usescratch=False)


#--- Long time-scale (complex) gain calibration for primary and secondary
#    phase (P) calibrators - 1037-295 and Titan: make table G2
#    Use the B-calibration on the fly (use 'gaintable' parameter).
gaincal(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.G2', spw =
	'*:16~112', field = '1037*,Titan', minsnr=1.0,
	solint= 'inf', selectdata=True, solnorm=False, refant = 'DV04',
	gaintable = 'cal-ngc3256.B1', calmode = 'ap')



#--- Inspect the G2 calibration table

#=== Commented out in production-run script
'''

plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'phase',
	poln='X', plotsymbol='o', plotrange = [0,0,-180,180], iteration= 'spw', 
        figfile='cal-phase_vs_time_XX.G2.png', subplot = 221, showgui=True)

plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'phase',
	poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], iteration= 'spw', 
        figfile='cal-phase_vs_time_YY.G2.png', subplot = 221, showgui=True)

plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'amp',
	poln='X', plotsymbol='o', plotrange = [], iteration = 'spw',
	figfile='cal-amp_vs_time_XX.G2.png', subplot = 221, showgui=True)

plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'amp',
	poln='Y', plotsymbol='o', plotrange = [], iteration = 'spw',
	figfile='cal-amp_vs_time_YY.G2.png', subplot = 221, showgui=True)

'''


#--- Bootstrapping of the absolute flux density for the phase calibrator 
#    1037-295 from Titan 
fluxscale( vis="ngc3256_line.ms", caltable="cal-ngc3256.G2",
        fluxtable="cal-ngc3256.G2.flux", reference="Titan",
        transfer="1037*", refspwmap=[0,1,1,1])


#--- Apply the calibrations to the science target and the phase calibrator

applycal( vis='ngc3256_line.ms', flagbackup=False, field='NGC*,1037*',
	interp=['nearest','nearest'], gainfield = ['1037*', '1037*'],
	gaintable=['cal-ngc3256.G2.flux', 'cal-ngc3256.B1'])


#--- Inspect calibrated data: Show the spectrum of science target


#=== Commented out in production-run script

'''

plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='amp',
	ydatacolumn='corrected', selectdata=True, field='NGC*',
	averagedata=True, avgchannel='', avgtime='10000',
	avgscan=True, avgbaseline=True, coloraxis='spw', showgui=True)

'''

#--- Inspect primary P calibrator


#=== Commented out in production-run script


#    amplitude vs time

'''
plotms(vis='ngc3256_line.ms', xaxis='time', yaxis='amp',
	ydatacolumn='corrected', selectdata=True, field='1037*',
	averagedata=True, avgchannel='128', avgtime='',
	avgscan=False, avgbaseline=False, coloraxis='spw', showgui=True)
'''


#    phase vs time

'''
plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='phase',
	ydatacolumn='corrected', selectdata=True, field='1037*',
	averagedata=True, avgchannel='', avgtime='1e8',
	avgscan=True, avgbaseline=False, coloraxis='baseline', showgui=True)
'''

#++++++++++ The above inspection of P calibrator would find the corrupted data
# unless we flagged them in advance (what we did). 
# The following sequence - repetition of calibration steps with late flagged
# data is thus unnecessary now and is commented out:

'''

#--- Flag corrupted data

flagdata(vis='ngc3256_line.ms', mode='manual',
	timerange='2011/04/16/04:13:35~04:13:45', flagbackup=True)

# copy-paste last args of flagdata() call, we will change just some of them
tget(flagdata)

timerange='2011/04/16/05:21:13~05:21:19'
flagdata()
timerange='2011/04/16/04:16:40~04:16:49'
flagdata()
timerange='2011/04/16/04:14:00~04:17:10'; antenna='PM03'
flagdata()
timerange='2011/04/17/01:52:20~01:53:10'; antenna='DV10'
flagdata()
timerange='2011/04/17/00:35:30~01:20:20'; antenna='DV04'; spw='3'
flagdata()



#--- Re-initialise the "MODEL" column for calibrators after flagging

delmod('ngc3256_line.ms')


#--- Repeat all the calibration with data shrinked by flagging
gaincal(vis='ngc3256_line.ms', caltable='cal-ngc3256.G1n', 
        spw='*:40~80', field='1037*',
        selectdata=True, solint='int', refant='DV04', calmode='p')


bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1n', 
         gaintable = 'cal-ngc3256.G1n', timerange='<2011/04/16/15:00:00',
         field = '1037*', minblperant=3, minsnr=2, solint='inf', 
         combine='scan,obs',
         bandtype='B', fillgaps=1, refant = 'DV04', solnorm=True)

bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1n', 
         gaintable = 'cal-ngc3256.G1n', timerange='>2011/04/16/15:00:00',
         field = '1037*', minblperant=3, minsnr=2, solint='inf', 
         combine='scan,obs', 
         bandtype='B', fillgaps=1, refant = 'DV04', solnorm=True, append=True)

setjy(vis='ngc3256_line.ms', field='Titan', 
      standard='Butler-JPL-Horizons 2012', 
      spw='0,1,2,3',usescratch=False)

gaincal(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.G2n', spw =
	'*:16~112', field = '1037*,Titan', minsnr=1.0,
	solint= 'inf', selectdata=True, solnorm=False, refant = 'DV04',
	gaintable = 'cal-ngc3256.B1n', calmode = 'ap')

fluxscale(vis="ngc3256_line.ms", caltable="cal-ngc3256.G2n",
          fluxtable="cal-ngc3256.G2n.flux", reference="Titan",
          transfer="1037*", refspwmap=[0,1,1,1])


applycal(vis='ngc3256_line.ms', flagbackup=False, field='NGC*,1037*',
         interp=['nearest','nearest'], gainfield = ['1037*', '1037*'],
         gaintable=['cal-ngc3256.G2n.flux', 'cal-ngc3256.B1n'])





#--- Inspect the phase calibrator again
plotms(vis='ngc3256_line.ms', xaxis='time', yaxis='amp',
	ydatacolumn='corrected', selectdata=True, field='1037*',
	averagedata=True, avgchannel='128', avgtime='',
	avgscan=False, avgbaseline=False, coloraxis='spw', showgui=True)


plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='phase',
	ydatacolumn='corrected', selectdata=True, field='1037*',
	averagedata=True, avgchannel='', avgtime='1e8',
	avgscan=True, avgbaseline=False, coloraxis='baseline', showgui=True)

# now OK

#--- /end of repeted steps

'''


#--- Imaging of calibrators
os.system('rm -rf result-phasecal_cont*')

# Primary P-calibrator - Inverse FFT and deconvolution in-one
clean(vis='ngc3256_line.ms', imagename='result-phasecal_cont', field='1037*',
      spw='*:20~120', selectdata=True, mode='mfs', niter=500,
      gain=0.1, threshold='0.75mJy', psfmode='hogbom',
      interactive=False, mask=[62, 62, 67, 67], imsize=128,
      cell='1arcsec', weighting='briggs', robust=0.0, nterms=2)


#--- make image statistics
calstat=imstat(imagename='result-phasecal_cont.image.tt0', region='', 
               box='85,8,120,120')
rms=(calstat['rms'][0])
print '>> rms in phase calibrator image: '+str(rms)
calstat=imstat(imagename='result-phasecal_cont.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)


#--- View the primary P-calibrator: set output to PNG file
imview(raster={'file': 'result-phasecal_cont.image.tt0', 'colorwedge':T,
               'range':[-0.004, 0.250], 'scaling':-1.5, 'colormap':'Rainbow 2'},
       out='result-phasecal_map.png', zoom=1)


#--- Apply calibrations for Titan
applycal(vis='ngc3256_line.ms', flagbackup=False, field='Titan',
	interp=['nearest', 'nearest'], gainfield = ['Titan', '1037*'],
	gaintable=['cal-ngc3256.G2.flux', 'cal-ngc3256.B1'])


os.system('rm -rf result-ampcal_cont*')


#--- Titan - Inverse FFT and deconvolution in-one
clean(vis='ngc3256_line.ms', imagename='result-ampcal_cont', 
      field='Titan', spw='0:20~120,1:20~120', mode='mfs', niter=200, 
      threshold='5mJy', psfmode='hogbom', mask=[62, 62, 67, 67], imsize=128,
      cell='1arcsec', weighting='briggs', robust=0.0, interactive=False)


#--- make image statistics
calstat=imstat(imagename="result-ampcal_cont.image",
               region="",box="85,8,120,120")
rms=(calstat['rms'][0])
print ">> rms in amp calibrator image: "+str(rms)
calstat=imstat(imagename="result-ampcal_cont.image",region="")
peak=(calstat['max'][0])
print ">> Peak in amp calibrator image: "+str(peak)
print ">> Dynamic range in amp calibrator image: "+str(peak/rms)


#--- View Titan: set output to PNG file

imview(raster={'file': 'result-ampcal_cont.image', 'colorwedge':T,
               'range':[-0.02, 0.250], 'scaling':-1.5, 'colormap':'Rainbow 2'},
       out='result-ampcal_map.png', zoom=1)
os.system('rm -rf ngc3256_line_target.ms*')


#--- Extract just the science target  NGC 3256 into separate MS

split(vis='ngc3256_line.ms', outputvis='ngc3256_line_target.ms',
	field='NGC*')