# Incremental calibration demonstration script -- gmoellen (05Jan25)
#
# This script demos and tests incremental G calibration.  The data 
#  from the opacity/gain-curve and interpolation demos has been 
#  adopted for this demo.
#
# The data used in this script is simulated: 6 antennas, 8 hours, 1 spw,
#  1 channel, 60s integrations, full pol, 1 point source target, 2 point
#  source calibrators.  The reference frequency is 22.235 GHz.
#  Both calibrators are observed for 2 minutes each
#  between 13 minute target source scans.  The primary calibrator has a
#  flux density of 2.345 Jy, and is at RA=19deg,Dec=75deg; the
#  secondary calibrator is 0.357 Jy, and is at RA=181deg,Dec=75deg.  The
#  target is a 13 mJy point source 1 degree west of the secondary calibrator.
#  The sources are circumpolar (always up), but the primary is in a
#  different direction (and so has different opacity effects).
#  The only systematic errors introduced in the simulation are opacity
#  (0.2n at zenith), the gain curve, and continuously varying G.
#  Noise (rms=10mJy each integration) has also been added.
#
# The calibration demo consists of solving first for G on the strongest
#  calibrator, field 1.  These solutions are interpolated onto a cumulative
#  table for all sources, which is then applied.  The resulting (point-source)
#  flux densities are then reported.  Then, an incremental G solution
#  is solved for on field 3.  This solution is accumulated onto the existing
#  cumulative table (field 3 only), and the field 3 data again corrected,
#  and flux densities reported.  The resulting flux density for field
#  three more closely matches the correct value after the incremental
#  calibration because its true flux density is enforced by the second
#  G solve.  This is true even though the SNR of the incremental (and thus
#  cumulative) solution is lower due to the lower flux density of field 3.


include 'logger.g'
dl.purge(0);

msname:='incrdemo.ms';

# Fill data from UVFITS (otf_interp.uvfits)
include 'ms.g'
myms:=fitstoms(msfile=msname,fitsfile='incremental.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(filename=msname);
mycal.initcalset();

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

# MS tool will be used to calculate flux densities
myms:=ms(filename=msname);  

# Report pre-cal (incoherent) apparent flux densities
mfd0:=[myms.ptsrc(1)[1],myms.ptsrc(2)[1],myms.ptsrc(3)[1]]
print 'Pre-cal f.d.    = ', mfd0


# Initial G calibration using field 1
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1]');
mycal.setapply(type='TOPAC',opacity=0.2);
mycal.setapply(type='GAINCURVE');
mycal.setsolve(type='G',table='cal.gcal1',t=15*60);
mycal.solve();
mycal.plotcal('PHASE','cal.gcal1');

# Accumulate these solution onto all fields on 120s timescale
mycal.accumulate(tablein='',incrtable='cal.gcal1',tableout='cal.cG1',
		 field=[1,2,3],calfield=1,interp='linear',t=120)
mycal.plotcal('PHASE','cal.cG1');

# Correct all data with these solutions
mycal.reset();
mycal.setdata();
mycal.setapply(type='TOPAC',opacity=0.2);
mycal.setapply(type='GAINCURVE');
mycal.setapply(type='G',table='cal.cG1',interp='linear');
mycal.correct();

# Report flux densities with calibration so far
mfd1:=[myms.ptsrc(1)[1],myms.ptsrc(2)[1],myms.ptsrc(3)[1]]
print 'Init-cal f.d.   = ', mfd1


# solve for scan-based (t=0) incremental solution on field 3
mycal.reset()
mycal.setdata(msselect='FIELD_ID==3');
mycal.setapply(type='TOPAC',opacity=0.2);
mycal.setapply(type='GAINCURVE');
mycal.setapply(type='G',table='cal.cG1',interp='linear');
mycal.setsolve(type='G',table='cal.gcal2',t=0);
mycal.solve();
mycal.plotcal('PHASE','cal.gcal2');

# accumulate this new solution onto field 3 only, make new cumulative table
mycal.accumulate(tablein='cal.cG1',incrtable='cal.gcal2',tableout='cal.cG2',
		 field=3,calfield=3,interp='linear');
mycal.plotcal('PHASE','cal.cG2');

# Correct field 3 only with this new table
mycal.reset();
mycal.setdata(msselect='FIELD_ID==3');
mycal.setapply(type='TOPAC',opacity=0.2);
mycal.setapply(type='GAINCURVE');
mycal.setapply(type='G',table='cal.cG2',interp='linear',select='FIELD_ID==3');
mycal.correct();
mycal.done();

# Report final flux densities
mfd2:=[myms.ptsrc(1)[1],myms.ptsrc(2)[1],myms.ptsrc(3)[1]];
print 'Final-cal f.d.  = ', mfd2;

myms.done();