###
### This is the script obtained from the archive.
##

import re

if re.search('^3.4', casadef.casa_version) == None:
 sys.exit('ERROR: PLEASE USE THE SAME VERSION OF CASA THAT YOU USED FOR GENERATING THE SCRIPT: 3.4')


print "# Running clean."

# You have not specified a source Id, I will assume you want to clean the science target(s).

clean(vis = 'calibrated.ms',
  imagename = 'calibrated.ms.image.continuum.source2',
  field = '2', # NGC 1614
  spw = '1,2,3',
  mode = 'mfs',
  interactive = T,
  imsize = [360, 360],
  cell = '0.05arcsec',
  phasecenter = 2,
  weighting = 'briggs',
  robust = 0.5)
exportfits(imagename='calibrated.ms.image.continuum.source2.image',fitsimage='source2.continuum.fits')

myimagebase = 'calibrated.ms.image.continuum.source2'
impbcor(imagename=myimagebase+'.image', pbimage=myimagebase+'.flux', outfile=myimagebase+'.image.pbcor', overwrite=True) # perform PBcorr

################################################################
# We modify it 
##############################################################

# 3 sources: J0423-013, Ceres, NGC1614
# NGC1614= field 2
# 4 spw
# lambda = 0.4 mm 
# FOV 7.6 arcsec

# Estimate sensitivity
# ToS = 50 min 
# Freq= 680GHz 
# Nant= 
# ---> sensitivity in 5.25GHz (3spw) = 0.4 mJy

# plotms: field=2, amp vs freq, averaging in time
# identify that the line is in spw=0
# plotms: field=2, spw=0, amp vs uvwave, average chan, time, 
# identify the longest baseline = 900 klambda
# corresponding to the ~resolution 0.23
# ---> cellsize ~0.05 
# ---> imsize 2*7.6/0.05 = 304 the one used in QA2 is larger so it is OK!




#######################
# Continuum image
#######################

visname='NGC1614_Band9.ms'
listobs(visname)


# First clean interactive

os.system('rm -rf NGC1614_continuum*')
clean(vis = visname,
  imagename = 'NGC1614_continuum',
  field = '2', # NGC 1614
  spw = '1,2,3',
  mode = 'mfs',
  interactive = T,
  imsize = [360, 360],
  cell = '0.05arcsec',
  weighting = 'briggs',
  robust = 0.5,
  niter=100000)

# Define masks around the bright emission
# after the first iterations new bright components appear add boxes...
# stop after 400 iterations
# The residuals show still something...
# 8.4 e-4

# Clean with threshold (~ 4*0.5mJy)

os.system('rm -rf NGC1614_continuum_thresh*')
clean(vis = visname,
  imagename = 'NGC1614_continuum_thresh',
  field = '2', # NGC 1614
  spw = '1,2,3',
  mode = 'mfs',
  interactive = F,
  imsize = [360, 360],
  cell = '0.05arcsec',
  weighting = 'briggs',
  robust = 0.5,
  threshold='1mJy',
  niter=100000)
# draw one circular mask around the galax and go to the end of clean
# Total clean flux 0.195


# Clean with threshold changing weigthing and using the same mask as before
# and interactive False

os.system('rm -rf NGC1614_continuum_thresh_natural*')
clean(vis = visname,
  imagename = 'NGC1614_continuum_thresh_natural',
  field = '2', # NGC 1614
  spw = '1,2,3',
  mode = 'mfs',
  interactive = F,
  imsize = [360, 360],
  cell = '0.05arcsec',
  mask = 'NGC1614_continuum_thresh.mask',
  weighting = 'natural',
  robust = 0.5,
  threshold='1mJy',
  niter=100000)

# Total clean flux 0.220

os.system('rm -rf NGC1614_continuum_thresh_uniform*')
clean(vis = visname,
  imagename = 'NGC1614_continuum_thresh_uniform',
  field = '2', # NGC 1614
  spw = '1,2,3',
  mode = 'mfs',
  interactive = F,
  imsize = [360, 360],
  cell = '0.05arcsec',
  mask = 'NGC1614_continuum_thresh.mask',
  weighting = 'uniform',
  robust = 0.5,
  threshold='1mJy',
  niter=100000)


# Total clean flux 0.2449

### Compare results!!!
### opening the three images (note different beam size)
### natural: 0.29x0.23
### uniform: 0.244x0.17
### briggs: 0.26x0.20
### measure rms (check in the xterm the difference between the three)
### natural: 8e-4
### uniform: 3e-3
### briggs: 7.8 e-4


######
#### Channel image to be used for analysis.
#####

## The expected sensitivity in 10 km/s is 4 mJy

### Plotms to visualize just spw=0 xaxis=velocity
### no restfreq fixed, so the velocity of the galaxy appear to be zero
### transform add Frame=BARY and RestFreq=691473 (restfreq of the CO(5-6))
### change xaxis to channel to identify continuum channels in spw 0
### 10 to 40 

visname='NGC1614_Band9.ms'


# subtract continuum
# want_cont False otherwise it takes too much time!!!!

os.system('rm -rf NGC1614_Band9.ms.contsub')
uvcontsub(visname, field='2', fitspw='0:10~40,1,2,3')


# Clean spectral line.
# now the contsub dataset has only the target so field=''
# mode='velocity'
# outframe='BARY', restfreq='691.473GHz'
# width= 10 km/s (the natural spectral res is ~7km/s)
# To cover from 4400 to 5000 nchan=60
# remember the units in start and width
# interactive T to fix the mask
# threshold= 12mJy

os.system('rm -rf NGC1614_CO6-5.*')
clean(vis = visname+'.contsub',
  imagename = 'NGC1614_CO6-5',
  field = '', # NGC 1614
  spw = '0',
  mode = 'velocity',
  nchan =60,
  outframe='BARY',
  start='4400km/s',
  width='10km/s',
  interactive = T,
  threshold='12mJy',
  niter=10000,
  imsize = [360, 360],
  restfreq='691.473GHz',
  cell = '0.05arcsec',
  weighting = 'briggs',
  robust = 0.5)

# I would stop after 8 cycles
# rms=7mJy



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