#
# SPLINE/POLY solver demo and test script - gmoellen (03Nov04)
#
# A script to demo and test G spline and B poly solvers.  The data
#  used in this script is simulated: 6 antennas, 8 hours, 1 spw, 
#  63 channels, 60s integrations, RR-only, 1 point source target, 2 point 
#  source calibrators (all within a few degrees of each other).  
#  Both calibrators are observed for 2 minutes each
#  between 13 minute target source scans.   This is typical of IRAM
#  observing.  The first calibrator has a known flux density of
#  2.345 Jy; the second calibrator is 0.357 Jy, but is considered "unknown"
#  in the script.  The target is a 13 mJy point source. 
#  The only systematic errors introduced in the simulation are B and G, with
#  G continuously varying in time.  Noise has been added (10mJy per spw per
#  integration).

include 'logger.g'

# fill demo dataset from uvfits...
include 'ms.g'
msname:='splnpolydemo.ms'
myms:=fitstoms(msfile=msname,fitsfile='splnpolydemo.uvfits');
myms.done()

# True flux densities (for later comparison)
tfd:=[2.345, 0.013, 0.357];

# Launch calibrater tool, initialize scratch columns
include 'calibrater.g'
mycal:=calibrater(msname);
mycal.initcalset();

# Set flux density of primary calibrator
include 'imager.g'
myimgr:=imager(msname)
myimgr.setjy(fieldid=1,fluxdensity=2.345);

myms:=ms(msname);  # to calculate flux densities


# First solve using conventional solvers
#---------------------------------------
# solve for temporary G on strongest source
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1]',
              mode='channel',nchan=30,start=15,step=1);
mycal.setsolve(type='G',table='cal.g0',t=120,refant=1);
mycal.solve();

# solve for B, using temp G 
mycal.reset();
mycal.setdata(msselect='FIELD_ID==1')
mycal.setapply(type='G',table='cal.g0');
mycal.setsolve(type='B',table='cal.b',t=86400,refant=1);
mycal.solve();

# solve for full-up G using B 
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1,3]');
mycal.setapply(type='B',table='cal.b');
mycal.setsolve(type='G',table='cal.g1',t=120,refant=1);
mycal.solve();
mycal.plotcal('AMP','cal.g1');

mycal.fluxscale('cal.g1','cal.g2','Source1','Source3')
mycal.plotcal('AMP','cal.g2');

# apply calibration to everything
mycal.reset();
mycal.setdata();
mycal.setapply(type='B',table='cal.b');
mycal.setapply(type='G',table='cal.g2');
mycal.correct();

# Query ms for calibrated source flux densities
mfd:=[myms.ptsrc(1), myms.ptsrc(2), myms.ptsrc(3)]
mfd:=floor(100000*mfd+0.5)/100000;   # make fewer sig fig
efd:=mfd-tfd;
print "";
print "Conventional Solvers:"
print "--------------------"
for (i in 1:3) {
 print spaste('Field ',i,': True=',tfd[i],
                          ' Meas=',mfd[i],
                          ' Err=',efd[i],
                          ' (',floor(1000*abs(efd[i]/tfd[i])+0.5)/10,'%)');
}



# Now solve using spline and poly solvers
#----------------------------------------

# solve for G spline in phase (1200s timescale), from center channels
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1,3]',
              mode='channel',nchan=30,start=15,step=1);
mycal.setsolvegainspline(table='cal.gspline0',
                         splinetime=1200,mode='PHAS',refant=1,npointaver=2);
mycal.solve();

# solve for B poly, using G spline, first calibrator only
mycal.reset();
mycal.setdata(msselect='FIELD_ID==1')
mycal.setapply(type='GSPLINE',table='cal.gspline0');
mycal.setsolvebandpoly('cal.bpoly',bpnorm=T,
                       maskcenter=0,maskedge=0,
                       refant=1,degamp=15,degphase=15)
mycal.solve();

# solve for G spline phase (1200s timescale) using B poly and all channels
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1,3]');
mycal.setapply(type='BPOLY',table='cal.bpoly');
mycal.setsolvegainspline(table='cal.gspline1',
                         splinetime=1200,mode='PHAS',refant=1,npointaver=2);
mycal.solve();


# solve for T using G spline phase and B poly, to get flux density
#  of second calibrater
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1,3]');
mycal.setapply(type='BPOLY',table='cal.bpoly');
mycal.setapply(type='GSPLINE',table='cal.gspline1');
mycal.setsolve(type='T',table='cal.t',t=120,refant=1)
mycal.solve();
mycal.plotcal('AMP','cal.t');
mycal.fluxscale('cal.t','cal.t2','Source1','Source3')
mycal.plotcal('AMP','cal.t2');


# set fieldid=2 fluxdensity according to fluxscale result
#  (TBD: get fluxscale fluxdensity automatically)
myimgr.setjy(fieldid=3,fluxdensity=0.357);
myimgr.done();


# solve for G spline in amplitude, using B and G phase
#  (note use of preavg to reduce amplitude bias)
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1,3]');
mycal.setapply(type='BPOLY',table='cal.bpoly');
mycal.setapply(type='GSPLINE',table='cal.gspline1');
mycal.setsolvegainspline(table='cal.gspline1',
                         splinetime=10800,mode='AMP',
                         preavg=10800/8,refant=1);
mycal.solve();

# apply calibration to everything
mycal.reset();
mycal.setdata();
mycal.setapply(type='BPOLY',table='cal.bpoly');
mycal.setapply(type='GSPLINE',table='cal.gspline1');
mycal.correct();

mycal.done();

# Query ms for calibrated source flux densities
mfd:=[myms.ptsrc(1), myms.ptsrc(2), myms.ptsrc(3)]
mfd:=floor(100000*mfd+0.5)/100000;
efd:=mfd-tfd;
print "";
print "Spline/Poly Solvers:"
print "--------------------"
for (i in 1:3) {
 print spaste('Field ',i,': True=',tfd[i],
                          ' Meas=',mfd[i],
                          ' Err=',efd[i],
                          ' (',floor(1000*abs(efd[i]/tfd[i])+0.5)/10,'%)');
}


myms.done();