# Baseline-based calibration (M,MF) demonstration script
# gmoellen (05Aug10)


# Demo 1:  Basic usage of M to remove closure errors.
#  In this example, the continuum visibility function of 3C237 is 
#  treated as a closure error, and removed using M.  

# Cleanup any previous runs:
shell('rm -Rf 3c273.* ngc5921* 3c286.* cal.* test.*');

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

# Fill 3C273 data
include 'ms.g';
myms:=fitstoms('3c273.ms','/aips++/data/demo/3C273XC1.fits');
myms.summary(T);
myms.done();

include 'calibrater.g';
mycal:=calibrater('3c273.ms');

mycal.initcalset();                 # Initialize model to unpolarized
                                    #   unit point source

mycal.setdata();                    # select whole data set (all 3c273)

mycal.setsolve(type='M',            # Arrange to solve for M
	       t=0,                 #  (per scan)
	       table='test.M');
mycal.solve();

mycal.reset();
mycal.setapply(type='M',            # Arrange to apply M
	       table='test.M');
mycal.correct();

#include 'msplot.g';
#mymp:=msplot('3c273.ms');

#  Plots of amp & phase show visibility of unit point source in total
#   intensity (parallel hands).  Since the model visibilities (via
#   initcalset) are unpolarized, the cross-hand visibilities have not 
#   been changed.  

# Now specify a non-zero polarization

include 'imager.g';
myimgr:=imager('3c273.ms');          
myimgr.setjy(fieldid=1,fluxdensity=[1.0,0.1,0.1,0.0]);  # Q=U=0.1*I
myimgr.done();

mycal.reset();                      # Resolve w/ polarized model:

mycal.setsolve(type='M',            # Arrange to solve for M
	       t=0,                 #  (per scan)
	       table='3c273.Mpol');
mycal.solve();

mycal.reset();
mycal.setapply(type='M',            # Arrange to apply M
	       table='3c273.Mpol');
mycal.correct();

# At this point, the plots should show the modelled polarization, too.

# Finally, we can see that M solutions can remove the noise entirely;
#  the noise is a closure error, after all.  This is done by
#  solving M on the timescale of the data itself (i.e., per timestamp).

mycal.reset();                      # Resolve w/ polarized model:

mycal.setsolve(type='M',            # Arrange to solve for M
	       t=1,                 #  (per timestamp
	       table='3c273.Mnoise');
mycal.solve();

mycal.reset();
mycal.setapply(type='M',            # Arrange to apply M
	       table='3c273.Mnoise');
mycal.correct();

mycal.done();

# The corrected data is now precisely identical to the model.  Thus 
#   we see that use of M can be dangerous!  If used unwisely, it can 
#   remove all information contained in the data itself, both the
#   source visibility and the noise, and replace it with 
#   user-supplied information from the model.  This is otherwise known
#   as falsifying your data, so be careful!

# Demo #2:  Use of M to access G-independent B solutions.
#  In this example, M is used to estimate the time-dependent product
#  of gain and visibility, enabling a G-independent B solution
#  to be obtained.  This is useful if a suitable source model isn't
#  available for the typical "pre-B" G solution.  Here, since
#  a point source model is adequate for G solutions on the 
#  B calibrator, we compare the B's derived from both methods 
#  to show they are statistically the same.  We also demonstrate
#  the use of 'MF', the freq-dep version of M, or, "closure error
#  bandpass".  Images of a target source using these different
#  calibrations are made for comparison.

myms:=fitstoms('ngc5921.ms','/aips++/data/demo/NGC5921.fits');
myms.summary(T);
myms.done();

# Flag ACs:
include 'autoflag.g';
myaf:=autoflag('ngc5921.ms');
myaf.setdata();
myaf.setselect(autocorr=T);
myaf.run();
myaf.done();

mycal:=calibrater('ngc5921.ms');

myimgr:=imager('ngc5921.ms');
myimgr.setjy(fieldid=1);


mycal.setdata(msselect='FIELD_ID==1');  # select 3C286

# First do it traditionally, G then B
mycal.setsolve(type='G',                # Arrange to solve for G
	       t=1,                     #   per timestamp
	       table='3c286.G');
mycal.solve();

mycal.reset();
mycal.setapply(type='G',                # Arrange to apply G
	       table='3c286.G');
mycal.setsolve(type='B',                # Arrange to solve for B
	       t=86000,                 #   one solution
	       table='3c286.B1');
mycal.solve();


# Now do it with M
mycal.reset();

mycal.setsolve(type='M',                # Arrange to solve for M
	       t=1,                     #   per timestamp
	       table='3c286.M');
mycal.solve();

mycal.reset();
mycal.setapply(type='M',                # Arrange to apply M
	       table='3c286.M');
mycal.setsolve(type='B',                # Arrange to solve for B
	       t=86000,                 #   one solution
	       table='3c286.B2');
mycal.solve();

# Compare results:  (mean ratio of bandpasses should be Amp~1, Ph~0)
b1:=table('3c286.B1');
g1:=b1.getcol('GAIN');
b2:=table('3c286.B2');
g2:=b2.getcol('GAIN');
r:=g1[1,1,1,,]/g2[1,1,1,,];
b1.done();
b2.done();

amean:=mean(abs(r));
aerr:=stddev(abs(r));
pmean:=mean(arg(r))*180/pi;
perr:=stddev(arg(r))*180/pi;

print 'Amp:   ', amean, '+/-', aerr;
print 'Phase: ', pmean, '+/-', perr, '(deg)';

include 'pgplotter.g';
mypg:=pgplotter();
mypg.subp(2,1);
mypg.hist(abs(r),amean-5*aerr,amean+5*aerr,100,0);
mypg.lab('Amp Ratio','N',' ');
mypg.hist(arg(r)*180/pi,pmean-5*perr,pmean+5*perr,100,0);
mypg.lab('Phase Diff (deg)','N',' ');


# Also get a MF (freq-dep M, or "baseline-based bandpass")
mycal.reset();
mycal.setapply(type='M',                # Arrange to apply M
	       table='3c286.M');
mycal.setsolve(type='MF',               # Arrange to solve for MF
	       t=86000,                 #   one solution
	       table='3c286.MF');
mycal.solve();

# Compare baseline-based solutions with antenna-based solutions

b3:=table('3c286.MF');
g3:=b3.getcol('GAIN');
b3.done();

# Form baseline 1-2 from g1 solution
g1_12:=g1[,,1,,1]*conj(g1[,,1,,2]);  

mypg.clear();

# Plot baseline 1-2 (last index=2) of baseline based solution
mypg.plotxy(seq(63),abs(g3[1,,2]),plotlines=F,ptsymbol=2,linecolor=2);
# Plot antenna-based 1-2 solution
mypg.plotxy(seq(63),abs(g1_12[1,1,]),newplot=F,linecolor=3);

# Thus we see that the baseline-based solution is a bit noiser than
#   the more sensitive version formed from the antenna-based solutions.   


# Finally, we compare antenna-based calibration with baseline-based
#  calibration on the target source image  

# Obtain time-dep gain solution for phase calibrater
mycal.reset();
mycal.setdata(msselect='FIELD_ID==2');

# solve for G, using B
mycal.setapply(type='B',
	       table='3c286.B1');

mycal.setsolve(type='G',
	       t=0,
	       table='cal.G');
mycal.solve();

# solve for M, using MF
mycal.reset();
mycal.setapply(type='MF',
	       table='3c286.MF');

mycal.setsolve(type='M',
	       t=0,
	       table='cal.M');
mycal.solve();


# Correct data in traditional way (G and B)...
mycal.setdata(msselect='FIELD_ID==3');
mycal.reset();
mycal.setapply(type='B',
	       table='3c286.B1');
mycal.setapply(type='G',
	       table='cal.G');
mycal.correct();


# ...and image:

myimgr.setdata(fieldid=3, mode='channel',  # Select data for field 3, spectral
               spwid=1, nchan=42, start=9, #  window 1 and all channels
               step=1);
myimgr.setimage(nx=256, ny=256, stokes='I',          # Imaging parameters
                cellx='10arcsec', celly='10arcsec',
                start=9, step=1, nchan=42,
                mode='channel',fieldid=3);
myimgr.weight(type='natural');                       # Natural weighting

myimgr.clean(algorithm='clark', niter=1000,          # Image and deconvolve
	     model='ngc5921_GB.model',
	     image='ngc5921_GB.image',
	     residual='ngc5921_GB.residual');


# Correct data with M and MF...
mycal.reset();
mycal.setapply(type='MF',
	       table='3c286.MF');
mycal.setapply(type='M',
	       table='cal.M');
mycal.correct();

# ... and image:
myimgr.clean(algorithm='clark', niter=1000,          # Image and deconvolve
	     model='ngc5921_MB.model',
	     image='ngc5921_MB.image',
	     residual='ngc5921_MB.residual');

mycal.done();
myimgr.done();

include 'viewer.g';
dv.gui();

#  The images produced in these two different ways are virtually identical.
#    If there were significant closure errors in this data, the baseline-based
#    solutions would compensate for them, and the resulting image would be
#    better.  On the other hand, the antenna-based solution is more sensitive
#    and thus should be favored if closure errors are not thought to be
#    a problem.  If the calibrator were a lot weaker, then the M
#    solutions would contribute more noise to the final image.