# Pipeline originally developed by Russ Taylor in 2011
# Modified by Ishwara Chandra in 2018 
# Major improvements in flagging and self-calibration 
# Kindly refer to the paper: Ishwara-Chandra et al 2020  
# (paper under review as on Feb 2020)
# NOT tested for uGMRT band-2 (150 MHz) 
# Queries: ishwar@ncra.tifr.res.in
#
# TEST.FITS is the output of gvfits, which coverts GMRT LTA format to multi-source FITS
# 
# Please CHANGE channels and source fields as per your data.
# Also change clip parameters if you have much stronger calibrator and/or target
#
importuvfits(fitsfile="TEST.FITS",vis="multi.ms",antnamescheme="old")
#
ms='multi.ms'
splitms1='target1.ms'
# If more than one target, split separately.
# splitms2='target2.ms'
# splitms3='target3.ms'
#
# The parameters below are typical for 200 MHz, 4096 channels at band-3 (300- 500 MHz)
# Please change as required for your data based on your bandwidth and no. of channels.
# In BAND-4, recommended channels corresponding to ~ 560 MHz to ~ 800 MHz. The sensitivity drops sharply after 800 MHz
# In any case DO NOT use beyond 820 MHz.
#
# CHANGE field numbers for calibrators and sources as needed
#
flagspw = ""             # 
gainspw = '0:401~3200'   # central ~ 75% good channel range for calibration
splitspw =  '0:101~4000' # split channel range
specave = 4              # number of channels to average; suggested post-average BW (approx)
                         # 1.5 MHz at band-5; 0.7 MHz at band-4; 0.4 MHz at band-3; 0.1 MHz at band-2
timeave = '0s'           # time averaging
gainspw2 = '0:201~500'   # central good channels after split for self-cal; here 780 channels after split
#
bpassfield     = '0'     # field number of the bandpass calibratorar, add phasecal if strong enough
fluxfield      = '0'     # field number of the primary flux calibrator
secondaryfield = '1'     # field number of the gain calibrator
gaincals       ='0,1'    # All calibrators
kcorrfield     = '0'     # field number of the antenna-based delay calibrator
target         = '2'     # If more than one target, use target='2,3,4...', 
                         # BUT modify the split part of the pipeline to split each target separately
                         # GO TO LINE 281 to edit the target field and change 'target' to '2', '3',...
refant='15,25'           # two options, first and second under flexi option.
#
# Change limits as required
clipfluxcal =[0.0,100.0] # atleast twice the expected flux; only to remove high points
clipphasecal =[0.0,50.0] # atleast twice the expected flux; only to remove high points
cliptarget =[0.0,30.0]   # atleast four times the expected flux; only to remove high points
clipresid=[0.0,10.0]     # 10 times the rms for single channel and single baseline
uvracal='>1.0klambda'    # uvrange for gain calibration in main calibration
uvrascal='>0.75klambda'   # uvrange for gain calibration in self calibration
#
# Filenames for initial round of calibration
kcorrfile0 = ms+'.kcal0'
bpassfile0 = ms+'.bcal0'
gainfilep0 =  ms+'.gcalp0'
gainfile0 =  ms+'.gcal0'
fluxfile0 =  ms+'.fluxscale0'
#
# Filenames for final round of calibration
kcorrfile = ms+'.kcal2'
bpassfile = ms+'.bcal2'
gainfilep2 =  ms+'.gcalp2'
gainfile =  ms+'.gcal2'
fluxfile =  ms+'.fluxscale2'
#
#
pcycles=4              # number of phase-only self-calibration
apcycles=4             # number of amplitude and phase self-calibration
                       # IMAGE MAY APPEAR TO HAVE WORSENED AFTER FIRST ROUND OF A&P SELF-CAL.
                       # WILL IMPROVE BY 3RD ROUND - CHECK AP3 and AP4 (IGNORE AP1 and AP2)
doflag=True            #
psolint=8.0            # this transaltes to 4, 2, 1 and 0.5 min for each phase-only selfcal loop
apsolint=16.0          # this transaltes to  8, 4, 2 and 1 min for each a&p  selfcal loop
startthreshold=0.05    # start threshold to stop flux (mJy, will reduce by count subsequently)
startniter=1500        # start iterations (will double in each phase-selfcal and 4 times in each A&P loop)
imagesize=[5880,5880]  # should cover alteast up to null at lower part of the band
cellsize='1.5arcsec'   # should be atleast 3 pixels in minor axis
wproj=-1               # w projection, default autocalculate
#
default(flagdata)
flagdata(vis=ms,mode="clip", spw=flagspw,field=fluxfield, clipminmax=clipfluxcal,
        datacolumn="DATA",clipoutside=True, clipzeros=True, extendpols=False, 
        action="apply",flagbackup=True, savepars=False, overwrite=True, writeflags=True)
flagdata(vis=ms,mode="clip", spw=flagspw,field=secondaryfield, clipminmax=clipphasecal,
        datacolumn="DATA",clipoutside=True, clipzeros=True, extendpols=False, 
        action="apply",flagbackup=True, savepars=False, overwrite=True, writeflags=True)
# After clip, now flag using 'tfcrop' option for flux and phase cal tight flagging
flagdata(vis=ms,mode="tfcrop", datacolumn="DATA", field=gaincals, ntime="scan",
        timecutoff=4.0, freqcutoff=3.0, timefit="line",freqfit="poly",flagdimension="freqtime", 
        extendflags=False, timedevscale=4.0,freqdevscale=4.0, extendpols=False,growaround=False,
        action="apply", flagbackup=True,overwrite=True, writeflags=True)
# Now extend the flags (80% more means full flag, change if required)
flagdata(vis=ms,mode="extend",spw=flagspw,field=gaincals,datacolumn="DATA",clipzeros=True,
         ntime="scan", extendflags=False, extendpols=True,growtime=80.0, growfreq=70.0,growaround=False,
         flagneartime=False, flagnearfreq=False, action="apply", flagbackup=True,overwrite=True, writeflags=True)
# Now flag for target - moderate flagging, more flagging in self-cal cycles
flagdata(vis=ms,mode="clip", spw=flagspw,field=target, clipminmax=cliptarget,
        datacolumn="DATA",clipoutside=True, clipzeros=True, extendpols=False, 
        action="apply",flagbackup=True, savepars=False, overwrite=True, writeflags=True)
flagdata(vis=ms,mode="tfcrop", datacolumn="DATA", field=target, ntime="scan",
        timecutoff=4.0, freqcutoff=3.0, timefit="line",freqfit="poly",flagdimension="freqtime", 
        extendflags=False, timedevscale=5.0,freqdevscale=5.0, extendpols=False,growaround=False,
        action="apply", flagbackup=True,overwrite=True, writeflags=True)
# Now extend the flags (80% more means full flag, change if required)
flagdata(vis=ms,mode="extend",spw=flagspw,field=target,datacolumn="DATA",clipzeros=True,
         ntime="scan", extendflags=False, extendpols=True,growtime=80.0, growfreq=80.0,growaround=False,
         flagneartime=False, flagnearfreq=False, action="apply", flagbackup=True,overwrite=True, writeflags=True)
# Now summary
flagdata(vis=ms,mode="summary",datacolumn="DATA", extendflags=True, 
         name=vis+'summary.split', action="apply", flagbackup=True,overwrite=True, writeflags=True)
#
print "Calibrating measurement set %s" % ms

print " starting setjy for flux calibrator"
setjy(vis=ms, field = fluxfield, spw = flagspw, scalebychan=True)

# Phase only cal added - suggested and tested by Silpa Sasikumar
#

gaincal(vis=ms,caltable=gainfilep0,field=bpassfield,spw=gainspw,intent="",
         selectdata=True,timerange="",uvrange="",antenna="",scan="",
         observation="",msselect="",solint="int",combine="",preavg=-1.0,
         refant=refant,refantmode="flex",minblperant=5,minsnr=1.0,solnorm=False,
         gaintype="G",smodel=[],calmode="p",append=False,splinetime=3600.0,
         npointaver=3,phasewrap=180.0,docallib=False,callib="",gaintable=[''],
         gainfield=[''],interp=[],spwmap=[],parang=True)


print " starting initial gaincal -> %s" % gainfile0
gaincal(vis=ms, caltable = kcorrfile0, field = kcorrfield, spw = flagspw, 
        refant = refant,  minblperant = 6, solnorm = True,  gaintype = 'K', 
        gaintable =[gainfilep0], solint = '10min', combine = 'scan', minsnr=1.0,
        parang = True, append = False)
  
print " starting bandpass -> %s" % bpassfile0
bandpass(vis=ms, caltable = bpassfile0, field = bpassfield, spw = flagspw, minsnr=1.0,
         refant = refant, minblperant = 6, solnorm = True,  solint = 'inf', 
         bandtype = 'B', fillgaps = 8, gaintable = [gainfilep0, kcorrfile0], gainfield = kcorrfield, 
         parang = True, append = False)
print " starting gaincal -> %s" % gainfile0
gaincal(vis=ms, caltable = gainfile0, field = gaincals, spw = gainspw, 
        refant = refant, solint = '1.0min', minblperant = 5, solnorm = False,  
        gaintype = 'G', combine = '', calmode = 'ap', minsnr=1.0, uvrange=uvracal,
        gaintable = [gainfilep0,kcorrfile0,bpassfile0], gainfield = [kcorrfield,bpassfield],
        append = False, parang = True)

print " starting fluxscale -> %s", fluxfile0 
fluxscale(vis=ms, caltable = gainfile0, reference = [fluxfield], 
          transfer = [secondaryfield], fluxtable = fluxfile0, 
          listfile = ms+'.fluxscale.txt0',
          append = False)               
#
print " applying calibrations: primary calibrator"
applycal(vis=ms, field = fluxfield, spw = flagspw, selectdata=False, calwt = False,
    gaintable = [gainfilep0,kcorrfile0,bpassfile0, fluxfile0],
    gainfield = [bpassfield,kcorrfield,bpassfield,fluxfield],
    parang = True)

print " applying calibrations: secondary calibrators"
applycal(vis=ms, field = secondaryfield, spw = flagspw, selectdata = False, calwt = False,
    gaintable = [gainfilep0,kcorrfile0, bpassfile0, fluxfile0],
    gainfield = [bpassfield,kcorrfield, bpassfield,secondaryfield],
    parang= True)

print " applying calibrations: target fields"
applycal(vis=ms, field = target, spw = flagspw, selectdata = False, calwt = False,
    gaintable = [gainfilep0,kcorrfile0, bpassfile0, fluxfile0],
    gainfield = [bpassfield,kcorrfield, bpassfield,secondaryfield],
    parang= True)
# Change clipmax as required
default(flagdata)
flagdata(vis=ms,mode="clip", spw=flagspw,field=fluxfield, clipminmax=clipfluxcal,
        datacolumn="corrected",clipoutside=True, clipzeros=True, extendpols=False, 
        action="apply",flagbackup=True, savepars=False, overwrite=True, writeflags=True)
flagdata(vis=ms,mode="clip", spw=flagspw,field=secondaryfield, clipminmax=clipphasecal,
        datacolumn="corrected",clipoutside=True, clipzeros=True, extendpols=False, 
        action="apply",flagbackup=True, savepars=False, overwrite=True, writeflags=True)
# After clip, now flag using 'tfcrop' option for flux and phase cal tight flagging
flagdata(vis=ms,mode="tfcrop", datacolumn="corrected", field=gaincals, ntime="scan",
        timecutoff=4.0, freqcutoff=3.0, timefit="line",freqfit="line",flagdimension="freqtime", 
        extendflags=False, timedevscale=4.0,freqdevscale=4.0, extendpols=False,growaround=False,
        action="apply", flagbackup=True,overwrite=True, writeflags=True)
# now flag using 'rflag' option  for flux and phase cal tight flagging
flagdata(vis=ms,mode="rflag",datacolumn="corrected",field=gaincals, timecutoff=3.0, 
        freqcutoff=3.0,timefit="poly",freqfit="line",flagdimension="freqtime", extendflags=False,
        timedevscale=4.0,freqdevscale=4.0,spectralmax=500.0,extendpols=False, growaround=False,
        flagneartime=False,flagnearfreq=False,action="apply",flagbackup=True,overwrite=True, writeflags=True)
# Now extend the flags (70% more means full flag, change if required)
flagdata(vis=ms,mode="extend",spw=flagspw,field=gaincals,datacolumn="corrected",clipzeros=True,
         ntime="scan", extendflags=False, extendpols=False,growtime=80.0, growfreq=80.0,growaround=False,
         flagneartime=False, flagnearfreq=False, action="apply", flagbackup=True,overwrite=True, writeflags=True)
# Now flag for target - moderate flagging, more flagging in self-cal cycles
flagdata(vis=ms,mode="clip", spw=flagspw,field=target, clipminmax=cliptarget,
        datacolumn="corrected",clipoutside=True, clipzeros=True, extendpols=False, 
        action="apply",flagbackup=True, savepars=False, overwrite=True, writeflags=True)
flagdata(vis=ms,mode="tfcrop", datacolumn="corrected", field=target, ntime="scan",
        timecutoff=4.0, freqcutoff=4.0, timefit="poly",freqfit="line",flagdimension="freqtime", 
        extendflags=False, timedevscale=5.0,freqdevscale=5.0, extendpols=False,growaround=False,
        action="apply", flagbackup=True,overwrite=True, writeflags=True)
# now flag using 'rflag' option 
flagdata(vis=ms,mode="rflag",datacolumn="corrected",field=target, timecutoff=4.0, 
        freqcutoff=4.0,timefit="line",freqfit="line",flagdimension="freqtime", extendflags=False,
        timedevscale=5.0,freqdevscale=4.0,spectralmax=500.0,extendpols=False, growaround=False,
        flagneartime=False,flagnearfreq=False,action="apply",flagbackup=True,overwrite=True, writeflags=True)
# Now summary
flagdata(vis=ms,mode="summary",datacolumn="corrected", extendflags=True, 
         name=vis+'summary.split', action="apply", flagbackup=True,overwrite=True, writeflags=True)
#
print " Deleting existing model column"
clearcal(ms)
print "Calibrating measurement set %s" % ms

print " starting setjy for flux calibrator"
setjy(vis=ms, field = fluxfield, spw = flagspw, scalebychan=True)

print " starting initial gaincal -> %s" % gainfile

gaincal(vis=ms,caltable=gainfilep2,field=bpassfield,spw=gainspw,intent="",
         selectdata=True,timerange="",uvrange="",antenna="",scan="",
         observation="",msselect="",solint="int",combine="",preavg=-1.0,
         refant=refant,refantmode="flex",minblperant=5,minsnr=1.0,solnorm=False,
         gaintype="G",smodel=[],calmode="p",append=False,splinetime=3600.0,
         npointaver=3,phasewrap=180.0,docallib=False,callib="",gaintable=[''],
         gainfield=[''],interp=[],spwmap=[],parang=True)

gaincal(vis=ms, caltable = kcorrfile, field = kcorrfield, spw = flagspw, 
        refant = refant,  minblperant = 6, solnorm = True,  gaintype = 'K', 
        gaintable =[gainfilep2], solint = '10min', combine = 'scan', minsnr=1.0,
        parang = True, append = False)
  
print " starting bandpass -> %s" % bpassfile
bandpass(vis=ms, caltable = bpassfile, field = bpassfield, spw = flagspw, minsnr=1.0,
         refant = refant, minblperant = 6, solnorm = True,  solint = 'inf', 
         bandtype = 'B', fillgaps = 8, gaintable = [gainfilep2, kcorrfile], gainfield = kcorrfield, 
         parang = True, append = False)
     
print " starting gaincal -> %s" % gainfile
gaincal(vis=ms, caltable = gainfile, field = gaincals, spw = gainspw, 
        refant = refant, solint = '1.0min', minblperant = 5, solnorm = False,  
        gaintype = 'G', combine = '', calmode = 'ap', minsnr=1.0, uvrange=uvracal,
        gaintable = [gainfilep2,kcorrfile,bpassfile], gainfield = [kcorrfield,bpassfield],
        append = False, parang = True)

print " starting fluxscale -> %s", fluxfile 
fluxscale(vis=ms, caltable = gainfile, reference = [fluxfield], 
          transfer = [secondaryfield], fluxtable = fluxfile, 
          listfile = ms+'.fluxscale.txt2',
          append = False)               
#
print " applying calibrations: primary calibrator"
applycal(vis=ms, field = fluxfield, spw = flagspw, selectdata=False, calwt = False,
    gaintable = [gainfilep2, kcorrfile,bpassfile, fluxfile],
    gainfield = [bpassfield,kcorrfield,bpassfield,fluxfield],
    parang = True)

print " applying calibrations: secondary calibrators"
applycal(vis=ms, field = secondaryfield, spw = flagspw, selectdata = False, calwt = False,
    gaintable = [gainfilep2, kcorrfile, bpassfile, fluxfile],
    gainfield = [bpassfield,kcorrfield, bpassfield,secondaryfield],
    parang= True)

print " applying calibrations: target fields"
applycal(vis=ms, field = target, spw = flagspw, selectdata = False, calwt = False,
    gaintable = [gainfilep2, kcorrfile, bpassfile, fluxfile],
    gainfield = [bpassfield,kcorrfield, bpassfield,secondaryfield],
    parang= True)
# Now flag for target - moderate flagging, just before split
flagdata(vis=ms,mode="clip", spw=flagspw,field=target, clipminmax=cliptarget,
        datacolumn="corrected",clipoutside=True, clipzeros=True, extendpols=False, 
        action="apply",flagbackup=True, savepars=False, overwrite=True, writeflags=True)
# now flag using 'rflag' option 
flagdata(vis=ms,mode="rflag",datacolumn="corrected",field=target, timecutoff=4.0, 
        freqcutoff=4.0,timefit="line",freqfit="line",flagdimension="freqtime", extendflags=False,
        timedevscale=5.0,freqdevscale=4.0,spectralmax=500.0,extendpols=False, growaround=False,
        flagneartime=False,flagnearfreq=False,action="apply",flagbackup=True,overwrite=True, writeflags=True)
# Now summary
flagdata(vis=ms,mode="summary",datacolumn="corrected", extendflags=True, 
         name=vis+'summary.split', action="apply", flagbackup=True,overwrite=True, writeflags=True)
#
split(vis=ms, outputvis = splitms1, datacolumn='corrected', 
          field = target, spw = splitspw, keepflags=False, width = specave)
# For more targets, add with exact targets as below - target1.ms, target2.ms,...
#split(vis=ms, outputvis = splitms2, datacolumn='corrected', 
#          field = '3', spw = splitspw, keepflags=False, width = specave)
#
print "\n Cleaning up. Starting imaging..."
#start self-calibration cycles    
ms=splitms1
count=1
scmode='p'
#
tclean(vis=ms, imagename=ms+'.'+scmode+str(count-1),imsize=imagesize,cell=cellsize, selectdata=True, datacolumn="corrected", 
       phasecenter="",stokes="I",projection="SIN", specmode="mfs",nchan=-1,gridder="widefield",wprojplanes=wproj,
       aterm=True,pblimit=-1, deconvolver="mtmfs",nterms=2,smallscalebias=0.6,restoration=True,pbcor=False,
       weighting="briggs",robust=0.0,uvtaper=[],niter=int(startniter),gain=0.1,
       threshold=str(startthreshold/(count))+'mJy',cyclefactor=1.1,
       minpsffraction=0.05,maxpsffraction=0.8,usemask="auto-multithresh",pbmask=0.0,sidelobethreshold=2.0,
       growiterations=75,restart=True,savemodel="modelcolumn",calcres=True,calcpsf=True,parallel=False)
exportfits(imagename=ms+'.'+scmode+str(count-1)+'.image.tt0', fitsimage=ms+'.'+scmode+str(count-1)+'.fits')
for j in range(pcycles):  
  scmode='p'
  pniter=int(startniter*(count+2)*2**count)
  if pniter >=160000:
     pniter=160000
  if(doflag==True and count>=1):
     flagdata(vis=ms,mode="clip", spw="",field='', clipminmax=clipresid,
              datacolumn="RESIDUAL",clipoutside=True, clipzeros=True, extendpols=False, 
              action="apply",flagbackup=True, savepars=False, overwrite=True, writeflags=True)
     flagdata(vis=ms,mode="rflag",datacolumn="RESIDUAL",field='', timecutoff=5.0, 
              freqcutoff=5.0,timefit="line",freqfit="line",flagdimension="freqtime", extendflags=False,
              timedevscale=5.0,freqdevscale=4.0,spectralmax=500.0,extendpols=False, growaround=False,
              flagneartime=False,flagnearfreq=False,action="apply",flagbackup=True,overwrite=True, writeflags=True)
     flagdata(vis=ms,mode="summary",datacolumn="RESIDUAL", extendflags=False, 
              name=ms+'temp.summary', action="apply", flagbackup=True,overwrite=True, writeflags=True)
#
  gaincal(vis=ms,caltable=ms+'.'+scmode+str(count),selectdata=False,solint=str(psolint/2**count)+'min',refant=refant,refantmode="flex",
        minblperant=6, spw=gainspw2,minsnr=1.0,solnorm=True,gaintype="G",calmode=scmode,append=False, uvrange=uvrascal, parang=True)
#
  applycal(vis=ms, selectdata=False,gaintable=ms+'.'+scmode+str(count), parang=True,calwt=False,applymode="calflag",flagbackup=True)  
#
  tclean(vis=ms, imagename=ms+'.'+scmode+str(count),imsize=imagesize,cell=cellsize, selectdata=True, datacolumn="corrected", 
       phasecenter="",stokes="I",projection="SIN", specmode="mfs",nchan=-1,gridder="widefield",wprojplanes=wproj,
       aterm=True,pblimit=-1, deconvolver="mtmfs",nterms=2,smallscalebias=0.6,restoration=True,pbcor=False,
       weighting="briggs",robust=0.0,uvtaper=[],niter=int(pniter),gain=0.1,
       threshold=str(startthreshold/(count))+'mJy',cyclefactor=1.1,
       minpsffraction=0.05,maxpsffraction=0.8,usemask="auto-multithresh",pbmask=0.0,sidelobethreshold=2.0,
       growiterations=75,restart=True,savemodel="modelcolumn",calcres=True,calcpsf=True,parallel=False)
  exportfits(imagename=ms+'.'+scmode+str(count)+'.image.tt0', fitsimage=ms+'.'+scmode+str(count)+'.fits')
  count = count + 1
#
casalog.post("Completed phase only self-calibration, going to A&P calibration Cycle")
count=1
for j in range(apcycles):  
  scmode = 'ap'
  if count>= 4:
   sfactor=4
  else:
   sfactor=count
  apniter=int(startniter*count*4*2**count)
  if apniter >=160000:
     apniter=160000
#
  if(doflag==True):
     flagdata(vis=ms,mode="clip", spw="",field='', clipminmax=clipresid,
              datacolumn="RESIDUAL",clipoutside=True, clipzeros=True, extendpols=False, 
              action="apply",flagbackup=True, savepars=False, overwrite=True, writeflags=True)
     flagdata(vis=ms,mode="rflag",datacolumn="RESIDUAL",field='', timecutoff=5.0, 
              freqcutoff=5.0,timefit="line",freqfit="line",flagdimension="freqtime", extendflags=False,
              timedevscale=4.0,freqdevscale=4.0,spectralmax=500.0,extendpols=False, growaround=False,
              flagneartime=False,flagnearfreq=False,action="apply",flagbackup=True,overwrite=True, writeflags=True)
     flagdata(vis=ms,mode="summary",datacolumn="RESIDUAL", extendflags=False, 
              name=ms+'temp.summary', action="apply", flagbackup=True,overwrite=True, writeflags=True)
#
  gaincal(vis=ms,caltable=ms+'.'+scmode+str(count),selectdata=False,solint=str(apsolint/2**sfactor)+'min',refant=refant,
        refantmode="flex",spw=gainspw2,minblperant=6, minsnr=1.0,solnorm=True,gaintype="G",calmode=scmode,append=False, parang=True)
#
  applycal(vis=ms, selectdata=False,gaintable=ms+'.'+scmode+str(count), parang=True,calwt=False,applymode="calflag",flagbackup=True)  
#
  tclean(vis=ms, imagename=ms+'.'+scmode+str(count),imsize=imagesize,cell=cellsize, selectdata=True, datacolumn="corrected", 
       phasecenter="",stokes="I",projection="SIN", specmode="mfs",nchan=-1,gridder="widefield",wprojplanes=wproj,
       aterm=True,pblimit=-1, deconvolver="mtmfs",nterms=2,smallscalebias=0.6,restoration=True,pbcor=False,
       weighting="briggs",robust=0.0,uvtaper=[],niter=apniter,gain=0.1,
       threshold=str(startthreshold/(2*count))+'mJy',cyclefactor=1.1,
       minpsffraction=0.05,maxpsffraction=0.8,usemask="auto-multithresh",pbmask=0.0,sidelobethreshold=2.0,
       growiterations=75,restart=True,savemodel="modelcolumn",calcres=True,calcpsf=True,parallel=False)
  exportfits(imagename=ms+'.'+scmode+str(count)+'.image.tt0', fitsimage=ms+'.'+scmode+str(count)+'.fits')
  count = count + 1
print "Completed processing AP \n"
#