# Imaging 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_Imaging_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/25, barta@asu.cas.cz, EU ARC, Czech node


#--- Show uv spectra for individual SPWs

# See the CO line in spw 0 and CN doublet with high- and low-frequency lines
# in spw 1 


#=== Commented out in production-run script

'''
plotms(vis='ngc3256_line_target.ms',spw='',xaxis='channel',yaxis='amp',
       avgtime='1e8',avgscan=True,iteraxis='spw', showgui=True)


'''
#--- SPW 2 has no line emission: Check here whether we can resolve
#    the source in continuum


#=== Commented out in production-run script

'''
plotms(vis='ngc3256_line_target.ms',spw='2',xaxis='uvdist',yaxis='amp',
       avgtime='1e8',avgchannel='128',coloraxis='baseline', showgui=True)
'''

#--- Do cleaning for CONTINUUM - select jsut the line-free frequency channels

# Remove prospective old version of NGC 3256 image
os.system('rm -rf result-ngc3256_cont*')


# Let CASA advise us on image scales
#=== Commented out in production-run script
'''

im.open('ngc3256_line_target.ms')
im.advise()
im.close()

# The 'im' object advises to use 3 arcsec pixelsize and 64x64 pixels image
# size. Do a bit oversampling and oversize the image with pixelsize=1arcsec
# and image size 100x100 pixels


# Deconvolution and clean: How do you proceed in reality - interactive with
# fitting the cleanmask based on the dirty-image emission.
# The 'spw' parameter selects just the line-free frequency channels - see the
# first plot above.

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=0, 
      imsize=100, 
      cell='1arcsec',
      weighting='briggs', 
      robust=0.0, 
      interactive=True, 
      usescratch=False)
'''

# Deconvolution and clean: Finalised version (based on interactive clean above)
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', mask=[42,43,59,60],
      imsize=100, cell='1arcsec', weighting='briggs', robust=0.0, 
      interactive=False, usescratch=False)


#--- Print image statistics

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


#--- 'Display' continuum image: Put the result into PNG file

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)


#--- Self-calibration 


# 'Flat' phase variations on an intermediate timescale (1/2 hour)
# Use just the continuum (line-free) freq. channels
gaincal(vis='ngc3256_line_target.ms', field='NGC*',
    	caltable='cal-ngc3256_cont_30m.Gp',
    	spw='0:20~53;71~120,1:70~120,2:20~120,3:20~120',
    	solint='1800s', refant='DV07', calmode='p',
    	minblperant=3)



# Show calibration tables: Phase variations with time for XX and YY polarisation

#=== Commented out in production-run script

'''
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, showgui=True)

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, showgui=True)

# Phase variations on intermediate time scale of 30 minutes are OK, let us
# continue with application of phase self-calibration.

'''


# Apply the self-calibration table

applycal(vis='ngc3256_line_target.ms', interp='linear',
         gaintable='cal-ngc3256_cont_30m.Gp')


#--- Re-clean the continuum after self-cal

os.system('rm -rf result-ngc3256_cont_sc1*')
clean( vis='ngc3256_line_target.ms', imagename='result-ngc3256_cont_sc1',
	spw='0:20~53;71~120,1:70~120,2:20~120,3:20~120', psfmode='hogbom', 
	mode='mfs', niter=500, threshold='0.13mJy', mask=[42,43,59,60], 
        imsize=100, cell='1arcsec', weighting='briggs', robust=0.0,
        interactive=False, usescratch=False)


#--- Make and print image statistics

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)


# Note the much lower noise in the image and correspondingly higher dynamic 
# range (10x !) compared with original image before self-cal.


#--- 'Display' re-calibrated continuum image: Put the result into PNG file

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



#--- Continuum subtraction: Use line-free channels for continuum fit
# Output MS written to ngc3256_line_target.ms.contsub
uvcontsub(vis = 'ngc3256_line_target.ms',
	  fitspw='0:20~53;71~120,1:70~120,2:20~120,3:20~120', solint ='int', 
	  fitorder = 1,
	  combine='spw')

#--- Clean line 

# Remove old versions of the image
os.system('rm -rf result-ngc3256_line_CO.*')

# A way how to find parameter 'restfreq' for the clean() task 
# (rest frequency for the CO line) directly in CASA. 
# Other way: E.g., open Splatalogue in ALMA OT and find the line there.

#=== Commented out in production-run script

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


'''

# Clean - interactive
# Reset the cleaning box (or better fitting mask) in order to fit emission
# for the respective channels. We can identify the box in which fits emission
# for all channels, or better, select the mask on per-channel basis. 
# Use iteration over channels in the viewer (the 'go one step forward' arrow).

#=== Commented out in production-run script

'''

clean(vis='ngc3256_line_target.ms.contsub', imagename='result-ngc3256_line_CO',
      spw='0:38~87', mode='channel', start='', nchan=-1, width='', 
      psfmode='hogbom', outframe='LSRK', restfreq='115.271201800GHz', 
      mask=[53,50,87,83], niter=500, interactive=True, imsize=128, 
      cell='1arcsec', 
      weighting='briggs', robust=0.0, threshold='5mJy',
      usescratch=False)

'''

# Clean - finalised for production run
# We use the large box that spans over emission in all channels here
# Better approach is to save cleaning mask per channel in iterative clean()
# above, kee p the file result-ngc3256_line_CO.mask, and use parameter 
# mask='result-ngc3256_line_CO.mask' in the subsequent non-interactive clean()
# call

clean(vis='ngc3256_line_target.ms.contsub', imagename='result-ngc3256_line_CO',
      spw='0:38~87', mode='channel', start='', nchan=-1, width='', 
      psfmode='hogbom', outframe='LSRK', restfreq='115.271201800GHz', 
      mask=[35,42,87,83], niter=500, interactive=False, imsize=128, 
      cell='1arcsec', 
      weighting='briggs', robust=0.0, threshold='5mJy',
      usescratch=False)


# Look at the image
# Select result-ngc3256_line_CO.image and type 'raster' in the dialog that opens

#=== Commented out in production-run script
'''
viewer()
'''