#
#
# On-the-fly interpolation demonstration script -- gmoellen (04Nov12)
#
# This script demos and tests ordinary G application with
#  interpolation.  The data from the opacity/gain-curve demo
#  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 
#  known 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,
#  but its flux density is considered "unknown" in the script.  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.
#

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

msname:='interpdemo.ms';

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

myms:=ms(filename=msname);  # Will be used to calculate final flux densities




# Now solve on fields 1 using opacity and gain curve:
#---------------------------------------
# solve for G
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.gcal',t=15*60);
mycal.solve();
mycal.plotcal('PHASE','cal.gcal');


# apply calibration to everything
mycal.reset();
mycal.setdata();
mycal.setapply(type='TOPAC',opacity=0.2);
mycal.setapply(type='GAINCURVE');
mycal.setapply(type='G',table='cal.gcal',interp='linear');
mycal.correct();

mycal.done();

# Query ms for calibrated source flux densities
mfd:=[myms.ptsrc(1)[1], myms.ptsrc(2)[1], myms.ptsrc(3)[1]]
mfd:=floor(10000*mfd+0.5)/10000;
efd:=mfd-tfd;
print "";
print "Using opac, gaincurve, and linear interp of G:" 
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();

# Retry with interp='nearest'

# Applying calibration with 'nearest' makes very little difference
#  for this dataset, as it turns out.  And for some solution intervals,
#  linear behaves somewhat worse than nearest, probably due to residuals 
#  evening out better (by accident) when using nearest.
#  Linear interpolation will be more important for
#  more widely separated calibration observations.  Note that increasing
#  the solution interval for this dataset is not a way to emulate this 
#  because the solutions so obtained will be biased due to data averaging
#  with coherence loss.