#  Read in the data files

include 'vlafiller.g';

vlafillerfromdisk(
    filename="AW602_B030502.xp1",
    msname="AW602.ms",
    overwrite=F, 
    project="AW602",
    bandname="Q",
    verbose=F,
    async=F);

vlafillerfromdisk(
    filename="AW602_B030502.xp2",
    msname="AW602.ms",
    overwrite=F, 
    project="AW602",
    bandname="Q",
    verbose=F,
    async=F);

vlafillerfromdisk(
    filename="AW602_B030502.xp3",
    msname="AW602.ms",
    overwrite=F, 
    project="AW602",
    bandname="Q",
    verbose=F,
    async=F);

vlafillerfromdisk(
    filename="AW602_C030508.xp1",
    msname="AW602.ms",
    overwrite=F, 
    project="AW602",
    bandname="Q",
    verbose=F,
    async=F);

vlafillerfromdisk(
    filename="AW602_C030508.xp2",
    msname="AW602.ms",
    overwrite=F, 
    project="AW602",
    bandname="Q",
    verbose=F,
    async=F);

#  Set the flux densities

include 'imager.g';
im:=imager('AW602.ms');

im.setjy(fieldid=-1, spwid=-1, fluxdensity=-1);
im.setdata(fieldid=11, spwid=-1);
im.close();
im.done();

#  Do some clipping and flagging

include 'flagger.g';
fl:= flagger (msfile='AW602.ms');
fl.reset();
fl.filter (
   operation='range',
   range=['0Jy','2Jy'],
   trial=F);

# Addition flagging needed

fl.reset();
fl.filter (
   operation='range',
   range=['0Jy','2Jy'],
   trial=F);

fl.reset()
fl.timerange (
      starttime="2003/05/02/21:27:10",
      endtime  ="2003/05/02/21:48:30",
      trial=F);

fl.reset()
   fl.setantennas (ants=10);
   fl.settimerange("2003/05/08/00:00:00","2003/05/10/00:00:00");
   fl.flag(F);

fl.reset()
   fl.setantennas (ants=24);
   fl.settimerange("2003/05/01/00:00:00","2003/05/10/00:00:00");
   fl.flag(F);

fl.reset()
   fl.setantennas (ants=19);
   fl.settimerange("2003/05/04/00:00:00","2003/05/10/00:00:00");
   fl.flag(F);

fl.reset()
   fl.setantennas (ants=22);
   fl.settimerange("2003/05/02/20:06:00","2003/05/02/20:35:00");
   fl.flag(F);

#  flag channels 1:3, 58:63

fl.reset()
   fl.setchan (1:3);
   fl.settimerange("2003/05/01/20:06:00","2003/05/10/20:35:00");
   fl.flag(F);

fl.reset()
   fl.setchan (58:63);
   fl.settimerange("2003/05/01/20:06:00","2003/05/10/20:35:00");
   fl.flag(F);

fl.done();

#  Do initial gain calibration

include 'calibrater.g';
c:=calibrater('AW602.ms');
c.reset();
c.setdata (
   mode = 'channel',
   nchan = 63,
   start = 1, step = 1,
   msselect = 'FIELD_ID IN [1,11,12]');
c.setsolve (
   type='G',
   t=60,
   append = F,
   refant=15,
   table='AW602_g1');
c.state();
c.solve();

#  Do the bandpass

c.reset();
c.setdata (
   mode = 'channel',
   nchan = 63,
   start = 1, step = 1,
   msselect = 'FIELD_ID==12 && SPECTRAL_WINDOW_ID IN [13,14,15,16]');
c.setapply (
   type='G',
   table='AW602_g1')
c.setsolve (
    type='B',
    t=10000000,
    refant=15,
    table='AW602_b1',
    append=F);
c.state();
c.solve();

#  Now do gain calibration with bandpass

c.reset();
c.setdata (
   mode = 'channel',
   nchan = 63,
   start = 1, step = 1,
   msselect = 'FIELD_ID IN [1,11,12]');

c.setapply (
   type='B',
   table='AW602_b1',
   spwmap=[13,14,15,16,13,14,15,16,13,14,15,16,13,14,15,16])
c.setsolve (
   type='G',
   t=20,
   append = F,
   refant=15,
   table='AW602_g2')
c.state();
c.solve();

#  Now apply the B and G calibrations

c.reset();
c.setdata (
   mode = 'channel',
   nchan = 63,
   start = 1, step = 1);

c.setapply (
   type='B',
   table='AW602_b1',
   spwmap=[13,14,15,16,13,14,15,16,13,14,15,16,13,14,15,16])
c.setapply (
   type ='G',
   table='AW602_g2',
   select='FIELD_ID==1',
   spwmap = [1,2,1,2,5,6,5,6,1,2,5,6,1,2,5,6]);
c.correct();

c.close();
c.done();

#  SHOULD PUT FLUXSCALE IN HERE, BUT I THINK I NEED
#  spwmap in this routine.


#  Final editing of strange large values

include 'flagger.g';
fl:= flagger (msfile='AW602.ms');
fl.reset();
fl.filter (
   operation='range',
   range=['0Jy','5Jy'],
   column='CORRECTED_DATA',
   trial=F);

fl.done();

#   DATA IS NOW EDITED AND CALIBRATED

#   MAKE A CLEAN IMAGE OF MAIN CALIBRATOR

include 'imager.g';
imgr:=imager('AW602.ms');

zname:="M0336";
zfieldid:=1;

imgr.setdata(fieldid=zfieldid, spwid=[1,2,5,6],
   nchan=63, start=1, step=1);

imgr.setimage(
   cellx='0.5arcsec', celly='0.5arcsec',
   nx=512, ny=512, stokes="I", spwid=[1,2,5,6], 
   fieldid=zfieldid);

imgr.weight(type='natural');

imgr.filter(
   type='gaussian',
   bmaj = '2.3arcsec',
   bmin = '2.3arcsec');

dmap:=spaste(zname, '.imap');
dbeam:=spaste(zname, '.ibem');
model:=spaste(zname, '.ipts');
clean:=spaste(zname, '.icln');
residual:=spaste(zname, '.ires');
complist:=spaste(zname, '.icom');

imgr.clean (
   algorithm="hogbom",
   niter=200, gain=0.1,
   displayprogress=F, 
   model=model,
   image=clean,
   residual=residual);

imgr.close();
imgr.done();

#  Trying to make dirty images of all target sources.
#  Channels 4 to 57, combining 3.

#  First the sources 2,3,4,5,6 at frequency 1

include 'imager.g';
imgr:=imager('AW602.ms');
for (n in 2:6)  {
  
zname:=spaste('FLD_F1_', n);
zfieldid:=n;
#  Set data
imgr.setdata(fieldid=zfieldid, spwid=[3,4],
   mode='channel',
   nchan=63, start=1, step=1);

#  Set image
imgr.setimage(
   cellx='0.5arcsec', celly='0.5arcsec',
   nx=512, ny=512, stokes="I", spwid=[3,4], 
   fieldid=zfieldid,
   nchan=54, start=4, step=3);

dmap:=spaste(zname, '.imap');
dbeam:=spaste(zname, '.ibem');


imgr.makeimage (
   type = 'corrected',
   image = dmap):
imgr.makeimage (
   type = 'psf',
   image = dbeam);

}

#  next the sources 2,7,8,9,10 at frequency 2

for (n in [2,7,8,9,10))  {
  
zname:=spaste('FLD_F2_', n);
zfieldid:=n;
#  Set data
imgr.setdata(fieldid=zfieldid, spwid=[7,8],
   mode='channel',
   nchan=63, start=1, step=1);

#  Set image
imgr.setimage(
   cellx='0.5arcsec', celly='0.5arcsec',
   nx=512, ny=512, stokes="I", spwid=[7,8], 
   fieldid=zfieldid,
   nchan=54, start=4, step=3);

dmap:=spaste(zname, '.imap');
dbeam:=spaste(zname, '.ibem');


imgr.makeimage (
   type = 'corrected',
   image = dmap):
imgr.makeimage (
   type = 'psf',
   image = dbeam);

}

imgr.close();
imgr.done();