listobs("sis14_twhya_calibrated_flagged.ms")

### comparison with without selfcal ###
os.system('rm -rf image_noselfcal.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms',
       imagename='image_noselfcal',
       field='5',
       spw='',
       specmode='mfs',
       deconvolver='hogbom',
       nterms=1,
       gridder='standard',
       imsize=[250,250],
       cell=['0.1arcsec'],
       weighting='natural',
       threshold='0mJy',
       niter=5000,
       interactive=False)

rms = 2.4859e-3
peak = 4.82125e-1


### continuum selfcal ###
os.system('rm -rf first_image.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms',
       imagename='first_image',
       field='5',
       spw='',
       specmode='mfs',
       deconvolver='hogbom',
       nterms=1,
       gridder='standard',
       imsize=[250,250],
       cell=['0.1arcsec'],
       weighting='natural',
       threshold='0mJy',
       niter=5000,
       interactive=True,
       savemodel='modelcolumn')

os.system("rm -rf phase.cal")
gaincal(vis="sis14_twhya_calibrated_flagged.ms",
        caltable="phase.cal",
        field="5",
        solint="30s",
        calmode="p",
        refant="DV22",
        gaintype="G")
applycal(vis="sis14_twhya_calibrated_flagged.ms",
         field="5",
         gaintable=["phase.cal"],
         interp="linear")

os.system("rm -rf sis14_twhya_selfcal.ms sis14_twhya_selfcal.ms.flagversions")
split(vis="sis14_twhya_calibrated_flagged.ms",
      outputvis="sis14_twhya_selfcal.ms",
      datacolumn="corrected")

os.system('rm -rf second_image.*')
tclean(vis='sis14_twhya_selfcal.ms',
       imagename='second_image',
       field='5',
       spw='',
       specmode='mfs',
       deconvolver='hogbom',
       nterms=1,
       gridder='standard',
       imsize=[250,250],
       cell=['0.1arcsec'],
       weighting='natural',
       threshold='0mJy',
       interactive=True,
       niter=5000,
       savemodel='modelcolumn')

os.system("rm -rf phase_2.cal")
gaincal(vis="sis14_twhya_selfcal.ms",
        caltable="phase_2.cal",
        field="5",
        solint="30s",
        calmode="p",
        refant="DV22",
        gaintype="G")

applycal(vis="sis14_twhya_selfcal.ms",
         field="5",
         gaintable=["phase_2.cal"],
         interp="linear")

os.system("rm -rf sis14_twhya_selfcal_2.ms sis14_twhya_selfcal_2.ms.flagversions")
split(vis="sis14_twhya_selfcal.ms",
      outputvis="sis14_twhya_selfcal_2.ms",
      datacolumn="corrected")

os.system('rm -rf third_image.*')
tclean(vis='sis14_twhya_selfcal_2.ms',
       imagename='third_image',
       field='5',
       spw='',
       specmode='mfs',
       deconvolver='hogbom',
       nterms=1,
       gridder='standard',
       imsize=[250,250],
       cell=['0.1arcsec'],
       weighting='natural',
       threshold='0mJy',
       interactive=True,
       niter=5000,
       savemodel='modelcolumn')

os.system("rm -rf amp.cal")
gaincal(vis="sis14_twhya_selfcal_2.ms",
        caltable="amp.cal",
        field="5",
        solint="30s",
        calmode="ap",
        refant="DV22",
        gaintype="G",
        solnorm=True)

applycal(vis="sis14_twhya_selfcal_2.ms",
         field="5",
         gaintable=["amp.cal"],
         interp="linear")

os.system("rm -rf sis14_twhya_selfcal_3.ms sis14_twhya_selfcal_3.ms.flagversions")
split(vis="sis14_twhya_selfcal_2.ms",
      outputvis="sis14_twhya_selfcal_3.ms",
      datacolumn="corrected")


os.system('rm -rf fourth_image.*')
tclean(vis='sis14_twhya_selfcal_3.ms',
       imagename='fourth_image',
       field='5',
       spw='',
       specmode='mfs',
       deconvolver='hogbom',
       nterms=1,
       gridder='standard',
       imsize=[250,250],
       cell=['0.1arcsec'],
       weighting='natural',
       threshold='0mJy',
       interactive=False,
       niter=5000)


### line imaging withouth selfcal ###
#os.system('rm -rf sis14_twhya_calibrated_flagged.ms.contsub')
uvcontsub(vis='sis14_twhya_calibrated_flagged.ms',fitspw='0:10~200;300~380', fitorder=0)

restfreq = '372.67249300GHz'
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
       imagename='n2hplus_noselfcal',
       field='5',
       spw='',
#       specmode='mfs',
       specmode='cube',
       start=230,width=1,nchan=20, 
       outframe='lsrk',veltype='radio', 
       restfreq=restfreq, 
       deconvolver='hogbom',
       nterms=1,
       gridder='standard',
       imsize=[250,250],
       cell=['0.1arcsec'],
       weighting='natural',
       threshold='0mJy',
       niter=5000,
       interactive=False)


### line imageing with continuum selfcal applied ###
spwmap_line = [0]
applycal(vis='sis14_twhya_calibrated_flagged.ms.contsub',
         spwmap=[spwmap_line, spwmap_line, spwmap_line], # entering the appropriate spwmap_line value for each spw in the input dataset
         field='5',
         gaintable=['phase.cal','phase_2.cal','amp.cal'],
         interp=['linear','linear', 'linear'])

split(vis='sis14_twhya_calibrated_flagged.ms.contsub',
      outputvis='sis14_twhya_calibrated_flagged.ms.contsub.selcal',
      datacolumn='corrected')


restfreq = '372.67249300GHz'
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub.selcal',
       imagename='n2hplus_selfcal',
       field='5',
       spw='',
#       specmode='mfs',
       specmode='cube',
       start=230,width=1,nchan=20, 
       outframe='lsrk',veltype='radio', 
       restfreq=restfreq, 
       deconvolver='hogbom',
       nterms=1,
       gridder='standard',
       imsize=[250,250],
       cell=['0.1arcsec'],
       weighting='natural',
       threshold='0mJy',
       niter=5000,
       interactive=False)