#
#
# OPACITY/GAINCURVE corrections demonstration script -- gmoellen (03Dec02)
#
# This script demos and tests ordinary G solving with and 
#  without opacity and gain curve corrections, illustrating the importance of
#  the these corrections in the determination of the flux density
#  scale.
# 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:='opacgcdemo.ms';

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


# First solve without using opacity / gaincurve
#----------------------------------------------
# solve for G
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1,3]');
mycal.setsolve(type='G',table='cal.g',t=60,refant=1);
mycal.solve();
mycal.plotcal('AMP','cal.g');

# Note that the amplitude solns contain variation due to opacity & gain curve.

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

# Note that the fluxscale result is consistent with the presumption
#  that the _average_ amplitude gain is the same for both calibrators, 
#  but that the difference varies substantially with time.  This cal
#  table will correct the gain variations well, but violates the fact
#  that the instrumental gains vary smoothly in time.  Flux densities
#  derived from this table will be incorrect.

# apply calibration to everything
mycal.reset();
mycal.setdata(msselect='FIELD_ID==1');
mycal.setapply(type='G',table='cal.g2',select='FIELD_ID==1');
mycal.correct();
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [2,3]');
mycal.setapply(type='G',table='cal.g2',select='FIELD_ID==3');
mycal.correct();

# 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 "Ignoring opacity/gaincurve corrections:"
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 opacity:
#------------------------
# solve for G
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1,3]');
mycal.setapply(type='TOPAC',t=-1,opacity=0.2);
mycal.setsolve(type='G',table='calopac.g',t=60,refant=1);
mycal.solve();
mycal.plotcal('AMP','calopac.g');

# Note that the amplitude solutions variations are considerably
#  less; the opacity effect has been accounted for by TOPAC.

mycal.fluxscale('calopac.g','calopac.g2','Source1','Source3')
mycal.plotcal('AMP','calopac.g2');

# Using the opacity corrections yields a fluxscaled cal table which
#  is much more consistent with the assumption that the antenna gains
#  are reasonably constant.  Note that the SNR is different for the 
#  two calibrators.  The G solution is still absorbing the
#  systematic gain curve errors (which are considerably smaller than
#  the opacity effects), and so the flux density scale is still not
#  quite right.  

# apply calibration to everything
mycal.reset();
mycal.setdata(msselect='FIELD_ID==1');
mycal.setapply(type='TOPAC',t=-1,opacity=0.2);
mycal.setapply(type='G',table='calopac.g2',select='FIELD_ID==1');
mycal.correct();
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [2,3]');
mycal.setapply(type='TOPAC',t=-1,opacity=0.2);
mycal.setapply(type='G',table='calopac.g2',select='FIELD_ID==3');
mycal.correct();

# 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 opacity (only) corrections:" 
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 opacity and gain curve:
#---------------------------------------
# solve for G
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [1,3]');
mycal.setapply(type='TOPAC',t=30,opacity=0.2);
mycal.setapply(type='GAINCURVE',t=30);
mycal.setsolve(type='G',table='calopacgc.g',t=60,refant=1);
mycal.solve();
mycal.plotcal('AMP','calopacgc.g');

# Now, the amplitude solutions do not contain any variation due 
#  to opacity or gain curve.

mycal.fluxscale('calopacgc.g','calopacgc.g2','Source1','Source3')
mycal.plotcal('AMP','calopacgc.g2');

# Now, the only systematic variations in the amplitude gains are 
#  the random-walk-like variations expected from the electronics
#  or variable troposphere.  Gross opacity effects and the gain curve
#  have been properly compensated, and the flux density scale has thus
#  been more accurately transferred to the secondary calibrator.
#  If the remaining amplitude variations are due to a variable troposphere,
#  the calibration transfer will be improved further by using a time-dependent
#  model for the opacity, e.g., based on system temperature and weather
#  data (TBD).

# apply calibration to everything
mycal.reset();
mycal.setdata(msselect='FIELD_ID==1');
mycal.setapply(type='TOPAC',t=30,opacity=0.2);
mycal.setapply(type='GAINCURVE',t=30);
mycal.setapply(type='G',table='calopacgc.g2',select='FIELD_ID==1');
mycal.correct();
mycal.reset();
mycal.setdata(msselect='FIELD_ID IN [2,3]');
mycal.setapply(type='TOPAC',t=30,opacity=0.2);
mycal.setapply(type='GAINCURVE',t=30);
mycal.setapply(type='G',table='calopacgc.g2',select='FIELD_ID==3');
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 opacity and gaincurve corrections:" 
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();