#
# Author: Steven T. Myers (NRAO) 2001-12-12
# Set of functions to apply opacity and gain correction to vla data
#
# Contents: 
#            gaincor(par)          apply gain correction to DATA column
#            opac(tau,el)          calculate opacity correction at vector el
#            gaincv(a1,a2,el)      gain curve for vecs an1,an2,el             
#
# Usage: 
#
# 1. define parameters in record par
# 2. if needed, make ms from fits using makemsf(par) from automap.g
# 3. run gaincor(par)
#

include 'measures.g';
include 'table.g';
include 'quanta.g';
#
# Example set up of paramters
#
par.msf := '20010204AK.ms';
par.tau0 := 0.10;
#
# Function to apply gain and opacity correction to ms
#
gaincor := function( par ) {
#
   msf := par.msf;
   tau0 := par.tau0;
# 
# Set the frame to be the appropriate telescope
#
   onam := spaste(msf, '/OBSERVATION');
   otab := table(onam);
   tel := otab.getcell("TELESCOPE_NAME",1);
   printf(' Data is from telescope %s\n',tel);
   telpos:=dm.observatory(tel);
   dm.measure(telpos,'WGS84');         # convert to geodetic (why????)
   dm.doframe(telpos);
   otab.close();
#
# Get source directions
#
   fnam := spaste(msf, '/FIELD');
   ftab := table(fnam);
   nfld := ftab.nrows();
   printf(' Found %d fields\n',nfld);
   fdirs := ftab.getcol("REFERENCE_DIR");
   fields := ftab.getcol("NAME");
   ftab.close();
#
# Process directions - store in record of directions
#
   dirq := [=];
   for ( i in 1:nfld ) {
      raq := dq.quantity(fdirs[1,,i],'rad');
      deq := dq.quantity(fdirs[2,,i],'rad');
      dirq[i] := dm.direction('j2000',raq,deq);
   }
#
# Make ms table
#
   mst := table(msf, readonly=F);
#
# Make ms table iterator and iterate
#
   ti := tableiterator(mst, "TIME");
   while ( ti.next() ) {
      st := ti.table();
      t := st.getcell("TIME",1);
      tq := dq.quantity(t,'s');
      tymd := dq.time(tq,form='ymd');
      print ' Processing timestamp ',tymd;
      tim := dm.epoch('tai',tymd) ; 
      dm.doframe(tim);
# proxy source as first row source for now, get EL for this set
      ifld := st.getcell("FIELD_ID",1)+1;   # (+1 corrects c++ indexing)
      azel := dm.measure(dirq[ifld],'AZEL');       
      el := azel.m1;                        # probably should use getvalue
#
# Get visibilties - this will be vector [shape=[4 1 nvis] ]
#
      nvis := st.nrows();
      vis := st.getcol("DATA");
#
# Get AN number vectors
#
      an1 := st.getcol("ANTENNA1")+1;
      an2 := st.getcol("ANTENNA2")+1;
#
# Construct el numeric vector
#
      elv := array( el.value, nvis );
#
# construct opacity corrections = ( 1 - tau/sin(EL) )^(-1)
#
      oc := opac( tau0, elv );
#
# construct gain corrections
#
      gc := gaincv( an1, an2, elv );
      vc := oc/gc;
#
# make compatible with vis vector
#
      vcv := array(1,4,nvis);
      for ( i in 1:nvis ) {
         vcv[1:4,i] := array( vc[i], 4 );
      }
# 
# apply to data vector
#
      vis *:= vcv;
#
# put back in ms
#
      st.putcol("DATA",vis);
      printf(' Corrected %s with %d rows at EL = %.2f\n',
	fields[ifld], nvis, dq.convert(el,'deg').value );
   }
#
# Done
#
   ti.done();
   mst.close();
   print 'Completed gain and opacity correction\n';
}
#
# Define the opacity correction
#
opac := function( tau0, el ) {
#
# Do atmospheric absorption correction on vector of el
# tau(el) = tau0/sin(el) -> oc = 1/[ 1 - tau(el) ]
# maximum tau = 1e3
   sel := sin( el );
   tau := tau0/sel;
   if ( tau > 1.e4 ) { tau := 1.e3 };
   oc := 1.0/( 1.0 - tau );
   return oc;
}
#
# Define a gain curve (would index these by antenna eventually)
#
gaincv := function( an1, an2, el ) {
# Example:
#  ANTEN 2 0; STOKES = 'L'; BIF=1; EIF=1; 
#  CLCOR 1.0041e+00, -6.0815e-04, 2.2647e-05, 0; 
# Could use polyfitter.eval function (in mathematics module) but
# is overkill for this
   c := [1.0041e+00, -6.0815e-04, 2.2647e-05 ];
   g := c[1] + c[2]*el + c[3]*el*el;
   if ( g < 1.e-4 ) { g := 1.e-4 };
   return g;
}