#############################################################################
#                     Script for the reduction of NGC 1333
##############################################################################
## Initialize
## Include relevant tools if not already done:
##
include 'vlafiller.g';
include 'ms.g';
include 'autoflag.g';
include 'autoflag.g';
include 'flagger.g';
include 'msplot.g';
include 'calibrater.g';
include 'viewer.g';
include 'imager.g';
include 'image.g';
include 'imagepol.g';
#
#
#############################################################################
## Create the MS file from the originbal FITS file
#############################################################################

mset:=fitstoms(msfile='ngc1333.ms',
       fitsfile='/home/rwarmels/aips++/tst1.1/wrkspace/data/N1333_1.UVFITS');
                                          #  in the AIPS++ Data Repository
mset.summary(verbose=T);                   # Obtain summary
mset.done();                               # Finish ms tool
#

##############################################################################
## Flag bad data
##############################################################################
#
af:=autoflag('ngc1333.ms');
af.setdata();
af.setselect(autocorr=T);
af.run();

## Flag the correlator gliches
af.setselect(spwid=1,field=1,timerng="2003/05/02/21:44 2003/05/02/21:50");
af.setselect(spwid=1,field=[2,3,4,5,6],timerng="2003/05/02/21:40 2003/05/02/21:51");
af.setselect(spwid=2,field=[7],timerng="2003/05/02/21:56 2003/05/02/21:58");
af.setselect(spwid=2,field=[8,9,10,11,12],timerng="2003/05/02/21:55 2003/05/02/22:20");
##
## Flag the non-fringing antennas: 9 and 14
af.setselect(ant=9,spwid=[1,2]);
af.setselect(ant=14,spwid=1,timerng="2003/05/02/18:36 2003/05/02/18:53");
af.setselect(ant=14,spwid=2,timerng="2003/05/02/18:52 2003/05/02/19:10");

##
## Flag Antenna 22, spwid=2
af.setselect(ant=22,spwid=2);
#
## Select the first and last three channels
af.setselect(chan=[1,3,61,63],spwid=[1,2]);
###
af.run();
af.done();
#

#########################################################################
# Set the flux densities
#########################################################################
imgr:=imager(filename='ngc1333.ms');

# All for Field Id = 1
imgr.setjy(fieldid=13);
imgr.setdata(fieldid=13,spwid=1);

# All for Field Id = 2
imgr.setjy(fieldid=14);
imgr.setdata(fieldid=14,spwid=2);

imgr.done()

#
##########################################################################
# Do the Gain, Flux and Band pass calibration
##########################################################################
## Derive gain calibration solutions
##
cal:=calibrater(filename='ngc1333.ms');
cal.reset();
cal.setdata(msselect='FIELD_ID IN [1,13,15] && SPECTRAL_WINDOW_ID IN [1]',
       mode='channel', start=5,  nchan=55);
cal.setapply(type='GAINCURVE', t=-1);
cal.setsolve(type='G', t=0, refant=26, table='ngc1333.1.gcal');
cal.solve();
cal.plotcal(tablename='ngc1333.1.gcal', multiplot=F);

cal.reset();
cal.setdata(msselect='FIELD_ID IN [7,14,16] && SPECTRAL_WINDOW_ID IN [2]',
           mode='channel', start=5,  nchan=55);
cal.setapply(type='GAINCURVE', t=-1);
cal.setsolve(type='G', t=0, refant=26, table='ngc1333.2.gcal');
cal.solve();
cal.plotcal(tablename='ngc1333.2.gcal', multiplot=F);

## Transfer the flux density scale:
cal.fluxscale(tablein='ngc1333.1.gcal', tableout='ngc1333.1.fluxcal',
       reference='0542+498_1',
       transfer=["0336+323_1 0319+415_1"]);
cal.plotcal(tablename='ngc1333.1.fluxcal', plottype='AMP', multiplot=F);

## Transfer flux scales
cal.fluxscale(tablein='ngc1333.2.gcal', tableout='ngc1333.2.fluxcal',
       reference='0542+498_2',
       transfer=["0336+323_2 0319+415_2"]);
cal.plotcal(tablename='ngc1333.2.fluxcal', plottype='AMP', multiplot=F);

## Derive bandpass calibration solutions:
cal.reset();
cal.setdata(msselect='FIELD_ID==15 && SPECTRAL_WINDOW_ID==1');
cal.setapply(type='GAINCURVE', t=-1);
cal.setapply(type='G', t=0.0, table='ngc1333.1.gcal');
cal.setsolve(type='B', t=86400.0, refant=26, table='ngc1333.1.bcal');
cal.solve();
cal.plotcal(tablename='ngc1333.1.bcal', plottype='AMP', multiplot=F);

cal.reset();
cal.setdata(msselect='FIELD_ID==16 && SPECTRAL_WINDOW_ID==2');
cal.setapply(type='GAINCURVE', t=-1);
cal.setapply(type='G', t=0.0, table='ngc1333.2.gcal');
cal.setsolve(type='B', t=86400.0, refant=26, table='ngc1333.2.bcal');
cal.solve();
cal.plotcal(tablename='ngc1333.2.bcal', plottype='AMP', multiplot=F);

## Correct gcal & source data
cal.reset();
cal.setdata(msselect='FIELD_ID IN [1,2,3,4,5,6] && SPECTRAL_WINDOW_ID==1')
cal.setapply(type='GAINCURVE', t=-1);
cal.setapply(type='G', t=0.0, table='ngc1333.1.fluxcal',
       select='FIELD_ID IN [1]');   # fld1 = gcal
cal.setapply(type='B', t=0.0,  table='ngc1333.1.bcal', select='');
cal.correct();

cal.reset();
cal.setdata(msselect='FIELD_ID IN [7,8,9,10,11,12] && SPECTRAL_WINDOW_ID==2')
cal.setapply(type='GAINCURVE', t=-1);
cal.setapply(type='G', t=0.0, table='ngc1333.2.fluxcal',
       select='FIELD_ID IN [7]');   # fld7 = gcal
cal.setapply(type='B', t=0.0,  table='ngc1333.2.bcal', select='');
cal.correct();

###########################################################################
# Done with the calibration
#
# split the data
#
include 'ms.g';
mset:= ms("ngc1333.ms",readonly=F);
mset.split('ngc1333.source.ms',fieldids=[2:6,8:12],spwids=[1,2],
whichcol='CORRECTED_DATA',nchan=[63,63],start=[1,1],step=[1,1]);
mset.done();
#
include 'autoflag.g';
af:=autoflag('ngc1333.source.ms');
af.setdata();
af.setselect(field=[2:6,8:12],clip=[expr="+ ABS RR LL",max=6]);
af.run();
af.done();
#
#
###########################################################################
# Now do the cleaning and get the independent channels
#
imgr:=imager('ngc1333.source.ms');
imgr.setdata(mode='channel',fieldid=[1:10],spwid=[1,2], nchan=[63,63],
   start=[1,1],step=[1,1]);
imgr.setimage(mode='channel',nx=640,ny=256,cellx='0.5arcsec',
   celly='0.5arcsec',stokes='I',nchan=18,start=3,step=5,fieldid=[2],
   spwid=[1,2]);
imgr.weight(type='briggs',rmode='norm',robust=0.5,mosaic=T);
imgr.setvp(dovp=T);
imgr.setoptions(ftmachine='mosaic',padding=1.0);
imgr.setmfcontrol(scaletype='SAULT',minpb=0.07,constpb=0.4)

imgr.makeimage(type='corrected',image='ngc1333.dirty.im')

# multi-scale CLEAN parameters
modname:='ngc1333.full.clean.md'
imgname:='ngc1333.full.clean.im'
resname:='ngc1333.full.clean.rs'
sclname:='ngc1333.full.clean.sc'

imgr.clean(algorithm='mfclark',niter=1000,
       gain=0.2,threshold='Jy',displayprogress=T,
   model=modname,image=imgname,residual=resname);
imgr.summary();
imgr.done()

im:=imager('ngc1333.full.clean.im');
im.moments(outfile='ngc1333.full.clean.mom0',moments=1,axis=4,mask='indexin(4,[4:12])', includepix=[0,100]);
im.moments(outfile='ngc1333.full.clean.mom1',moments=1,axis=4,mask='indexin(4,[4:12])', includepix=[0.009,100]);
im.done();