# Correction for source size? (3C48, generic if user supplied?).
  
# Need ability to make use of this in other calibrations.
# 1) replace filled Tcal values with calculated ones.
#    a) per receiver, all scans
#    a) per receiver, limited scans
#    c) per receiver, time range
# 2) calibrate with a given Tcal in the argument list, default to
#     filled Tcal.
# 3) add additional unit of just Ta so that no correction for eta
#     and tau are done there.
# 4) Direct use of Scal to get units of Jy.
# 5) Improved ability to specify tau and eta so that agreement between
#    Scal and calib can be assured.
# 6) Integrate these results seamlessly into Jim's K-band calibration scripts.

# priority 6, 4 & 2, 5, 3, 1

include 'ott.g';

# onscan = scan number for on-source data
# offscan = scan number for off-source data.  Uses LASTOFF from GO FITS file if not specified.
# fres = frequency resolution (MHz) to use in solving for Scal.  Defaults to 1 Mhz.
# calsource = source name of calibrator.  If unset use source name in onscan.
# calflux = flux (Jy) of calsource.  If unset, get from Ott paper (ott.g).
# eta = aperture efficiency.
# tauzenit = opacity at zenith
# plot = if T, will plot the scal_vec spectra
# overlay = if T and plot is T, will overlay plot of scal_vec on dish plotter's current plot.

# returns a record containing 
#   scal (average, one per polarization),
#   scal_vec (one per polarization, at specified resolution)
#   tcal (average, one per polarization),
#   tcal_vec (one per polarization, at specified resolution)
#   freq (frequency values used to calculate scal_vec and tcal_vec)
#   onscan - on source scan number, copy of input argument
#   offscan - actual off source scan number used
#   fres - actual frequency resolution used, in MHz
#   ns - number of channels that correspond to fres
#   elev - elevation (deg) used in atmosphere (tau) correction
#   calflux - the calibrator flux at each freq value
#   tauzenit - the actual opacity at zenith used
#   eta - the actual aperture efficiency used.
#   receiver - the name of the receiver
#   feed - the name of the feed
#   receptors - the name of the receptors, one for each polarization
#   polarization - the name of the polarization,  one for each receptor

# This will return a fail (is_fail(result) will be T) if something
# goes wrong - usually a defaulted argument isn't sufficient to
# determine what the true value to be used is (e.g. the calsource
# isn't in available by that name, the offscan or onscan can't be
# found, etc).

getscal := function  (onscan, offscan=unset, fres=unset, calsource=unset,
		      calflux=unset, eta=unset, tauzenit=unset,
		      nif=F,nfeed=1,  plot=F, overlay=F) 
{
    result := [=];
    # get the switching state information for onscan

    q := d.qsp(onscan);
    if (is_boolean(q)) {
	fail('Could not find requested on scan');
    }
    if (len(q.uphases) !=2) {
	fail('There are more than 2 switching states in the on scan, this only works for cal-on and cal-off');
    }
    # get the data for sig with cal and sig without cal
    vsigon := d.getr(onscan, 2, nif=nif, nfeed=nfeed);
    vsigoff := d.getr(onscan, 1, nif=nif, nfeed=nfeed);

    if (is_fail(vsigon) || is_fail(vsigoff) || is_boolean(vsigon) || is_boolean(vsigoff)) {
	fail('Could not retrieve raw data for on scan');
    }

    # verify that vsigon is cal on and vsigoff is cal off
    if (vsigon.other.state.CAL == vsigoff.other.state.CAL) {
	# oops, same CAL state, bail out
	fail('Both states in the on scan have the same CAL switching state.');
    }
       if (vsigon.other.state.CAL == 0) {
	# reverse them
	tmp := vsigon;
	vsigon := vsigoff;
	vsigoff := tmp;
    }

    if (is_unset(offscan)) {
	if (has_field(vsigon.other,'gbt_go') &&
	    has_field(vsigon.other.gbt_go,'LASTOFF')) {
	    offscan := vsigon.other.gbt_go.LASTOFF;
	} else 	if (has_field(vsigon.other,'nrao_gbt_glish') &&
		    has_field(vsigon.other.nrao_gbt_glish,'LASTOFF')) {
	    offscan := vsigon.other.nrao_gbt_glish.LASTOFF;
	}
        if (is_unset(offscan) || offscan==0) {
	    # can not proceed.  Bail out now.
	    fail('offscan can not be determined from information in onscan');
	}
    }

    # get the switching state information for offscan
    q := d.qsp(offscan);
    if (is_boolean(q)) {
	fail(paste('Could not find requested off scan:',offscan));
    }
    if (len(q.uphases) !=2) {
	fail('There are more than 2 switching states in the off scan, this only works for cal-on and cal-off');
    }
    # get the data for ref with cal and ref without cal
    vrefon := d.getr(offscan, 2, F, F, nif=nif, nfeed=nfeed);
    vrefoff := d.getr(offscan, 1, F, F, nif=nif, nfeed=nfeed);

    if (is_fail(vrefon) || is_fail(vrefoff) || is_boolean(vrefon) || is_boolean(vrefoff)) {
	fail('Could not retrieve raw data for off scan');
    }

    # verify that vrefon is cal on and vrefoff is cal off
    if (vrefon.other.state.CAL == vrefoff.other.state.CAL) {
	# oops, same CAL state, bail out
	fail('Both states in the off scan have the same CAL switching state.');
    }
    if (vrefon.other.state.CAL == 0) {
	# reverse them
	tmp := vrefon;
	vrefon := vrefoff;
	vrefoff := tmp;
    }

    result.onscan := onscan;
    result.offscan := offscan;

    # sanity check on conformance between ref and sig spectra
    # only check shape, not actual frequency and polarization values (trust that the user knows
    # what they are doing there)
    if (vsigon.data.arr::shape != vrefon.data.arr::shape) {
	fail('On and off scan numbers have different data shapes');
    }

    if (is_unset(fres)) fres := 1.0;
    # nearest channel to requested resolution
    ns := as_integer(fres*1.0e6/vsigon.data.desc.chan_width + 0.5);
    # if (ns > vsigon.data.arr::shape[2]) ns := vsigon.data.arr::shape[2];
    # if (ns < 1) ns := 1;
    result.fres := ns * vsigon.data.desc.chan_width / 1.0e6;
    result.ns := ns;

    vsig := vsigoff
    vsig.data.arr := (vsigoff.data.arr + vsigon.data.arr)/2.;
    vref:= vrefoff;
    vref.data.arr := (vrefoff.data.arr + vrefon.data.arr)/2.;

    vsig.data.arr := (vrefon.data.arr - vrefoff.data.arr) / (vsig.data.arr - vref.data.arr);
  
    d.uniput("globalscan1",vsig);

    # this does decimation
    d.boxcar(ns);

    globalscan1 := d.uniget('globalscan1');
    a := d.getvfarray();

    if (is_unset(calsource)) calsource := vrefon.header.source_name;
    if (!any(ott.sources()==calsource)) {
	# not found
	fail(paste('source name',calsource,
		   'not found in Ott sources, type ott.sources() for complete list'));
    }

    if (is_unset(calflux) || calflux <= 0.0) {
	calflux := ott.sjy(calsource, a.f/1e6);
	if (any(calflux<=0.0)) {
	    # some of calflux must have been out of range
	    if (all(calflux<=0.0)) {
		# all out of range
		fail('The observed frequencies are all out of range in the Ott table for this source');
	    } else {
		# assume constant at the edges
		note('Some of the observed frequency band is out of range in the Ott table for this source.');
		note('Assuming a constant flux from the nearest frequency in the range.');
		lowzero := T;
		for (i in 1:len(calflux)) {
		    if (calflux[i] > 0.0 && lowzero) {
			if (i > 1) {
			    for (j in 1:(i-1)) {
				calflux[j] := calflux[i];
			    }
			}
			lowzero := F;
		    } else if (calflux[i] == 0.0 && !lowzero) {
			if (i < len(calflux)) {
			    for (j in (i+1):len(calflux)) {
				calflux[j] := calflux[i];
			    }
			}
			break;
		    }
		}
	    }
	}
    } else {
	# make constant vector of len(a.f) - one for each frequency in decimated average
	calflux := array(calflux, len(a.f));
    }

    elevpi := dq.getvalue(dq.convert(dm.getvalue(globalscan1.header.azel)[2],'rad')); 
    result.elev := dq.getvalue(dq.convert(dm.getvalue(globalscan1.header.azel)[2],'deg'));

    eff := F;

   avgf := sum(a.f/len(a.f));
   eff := d.eff(avgf, result.elev);

    # scal_vec gets same shape as globalscan1.data.arr
    scal_vec := globalscan1.data.arr;
    tcal_vec := globalscan1.data.arr;
    scal := array(0.0, scal_vec::shape[1]);
    nchan := scal_vec::shape[2];

    result.calflux := calflux;
    result.freq := a.f;

 if (is_unset(tauzenit)) {
       for (i in 1:nchan) {
            fr := result.freq[i]/1e9
            tauzenit_vec[i] := 120.419173 + -11.2108453*fr + 
			       0.391668037*fr^2 + -0.00608775352*fr^3 + 
	                       3.55561155e-05*fr^4
# tauzenit calcs for Q-band Tcal calcs y=120.419173-11.2108453*fr+0.391668037*fr^2
#					 -0.00608775352*fr^3+3.55561155e-05*fr^4  
#					(07 Jan 2005 and 23 May 2005)
# tauzenit calcs for K-band Tcal calcs y=165943.778-65532.4584x+11449.4684*x^2-1161.48234*x^3
#                                        +75.3882436*x^4-3.24654712*x^5+0.0927548503*x^6
#                                        -0.0016952095*x^7+1.79825603e-05*x^8-8.43500173e-08*x^9
#                                        (17 May 2005)
# tauzenit calcs for X-band HIGH Tcal calcs y=0.00025x+0.0055  (04 Feb 2005 and 09 May 2005)
# tauzenit calcs for X-band LOW Tcal calcs y=0.001x+0.001075  (10 Nov 2004 and 09 May 2005)
# tauzenit=0.0063 for S-band Tcal calcs (09 Feb 2005 and 09 May 2005)
# tauzenit=0.0063 for L-band Tcal calcs (27 Apr 2005)
       }
    } else {
       for (i in 1:nchan) {
           tauzenit_vec[i] := tauzenit;
       }
    }

    result.tauzenit_vec := tauzenit_vec
    result.tauzenit := sum(tauzenit_vec)/nchan;
    atmos := exp(-tauzenit_vec / sin(elevpi));

    for (i in 1:len(scal)) {
       for (j in 1:nchan) {
	scal_vec[i,j] := calflux[j] * atmos[j] * globalscan1.data.arr[i,j];
          }
      scal[i] := sum(scal_vec[i,]) / nchan;
    }
    
    result.scal := scal;
    result.scal_vec := scal_vec;

# print eta, eff.eta
    if (is_unset(eta)) {
       for (i in 1:nchan) {
            invwl2 := (result.freq[i]/1000000000./30.)^2
            eta_vec[i] := (1.66152638/2.8) * exp(-0.054304736 * invwl2)
       }
# eta_vec[i] := (1.66152638/2.8) * exp(-0.054304736 * invwl2) for K-band Tcal calcs (18 May 2005)
# eta_vec[i] := (1.68849866/2.8) * exp(-0.093893884 * invwl2) for K-band Tcal calcs (14 Nov 2004)
# eta = 0.70 for S-band and X-band Tcal calcs (09 Feb 2005)
# eta = 0.70 for L-band Tcal calcs (27 Apr 2005)
    } else {
       for (i in 1:nchan) {
           eta_vec[i] := eta
       }
    }
    
    result.eta_vec := eta_vec
    result.eta := sum(eta_vec)/nchan

    tcal := scal * 2.8 * result.eta
    result.tcal := tcal;

    for (i in 1:len(scal)) {
       for (j in 1:nchan) {
           tcal_vec[i,j] := scal_vec[i,j] * 2.8 * eta_vec[j];
        }
    }

    result.tcal_vecl := tcal_vec;

    vtsys := vrefon;
    vtsys.data.arr := vrefoff.data.arr/(vrefon.data.arr - vrefoff.data.arr)  
    d.uniput("globalscan1",vtsys);
    d.boxcar(ns);
    vtsys := d.uniget('globalscan1');    
    vtsys.data.arr := vtsys.data.arr*tcal_vec

    result.tsys_vec := vtsys.data.arr

    if (plot) {
	s := globalscan1;
	s.data.arr := result.scal_vec;
	d.plotscan(s,overlay);
    }

    # and add in feed descriptions
    if (has_field(globalscan1.other.feed,'GBT_FEED_NAME')) {
	result.feed := globalscan1.other.feed.GBT_FEED_NAME;
    } else {
	result.feed := '';
    }
    if (has_field(globalscan1.other.feed,'GBT_RECEIVER')) {
	result.receiver := globalscan1.other.feed.GBT_RECEIVER;
    } else {
	result.receiver := '';
    }
    cp := globalscan1.other.polarization.CORR_PRODUCT;
    pt := globalscan1.other.feed.POLARIZATION_TYPE;
    result.polarization := array('',globalscan1.data.arr::shape[1]);
    if (has_field(globalscan1.other.feed,'GBT_RECEPTORS')) {
	rcptrs := globalscan1.other.feed.GBT_RECEPTORS;
    } else {
	rcptrs := array('',len(pt));
    }
    result.receptors := result.polarization;
    if (len(result.polarization) == cp::shape[2]) {
	for (i in 1:len(result.polarization)) {
	    # assumes auto-correlation cp[1,i] == cp[2,i]
	    thisRcptr := cp[1,i] + 1;
	    if (thisRcptr <= len(pt) && thisRcptr > 0) {
		result.polarization[i] := pt[thisRcptr];
		result.receptors[i] := rcptrs[thisRcptr];
	    } # else leave it as is
	}
    } else if (len(result.polarization == len(pt))) {
	# must be older data
	result.polarization := pt;
	# won't have receptor information
    } # else I give up

    return result;
}
  
caldump := function ( data_res, fid, bdrop = 0, edrop = 0 ) {
    start := bdrop+1
    end := len(data_res.freq)-edrop
    for (i in seq(start,end,1)) {
       fprintf(fid, '%-6d %-6d %s  %-8.3f %-8.3f %-12.6f %-12.6f %-8.3f %-8.3f %-8.3f %-8.3f %-8.3f %-8.3f %-8.3f \n',
                i, data_res.onscan, data_res.feed,
                data_res.elev, data_res.calflux[i], data_res.freq[i]/1000000,  data_res.tauzenit_vec[i], data_res.eta_vec[i],
                data_res.scal_vec[1,i], data_res.scal_vec[2,i], data_res.tcal_vecl[1,i], data_res.tcal_vecl[2,i], 
                data_res.tsys_vec[1,i], data_res.tsys_vec[2,i])
  }
}

findResids := function ( onscan ) {

     result := [=];

     vsigon := d.getr(onscan, phase=2, nfeed=1, nif=nif)
     vsigoff := d.getr(onscan, phase=1, nfeed=1, nif=nif)
     vrefon := d.getr(onscan-1, phase=2, nfeed=1, nif=nif)
     vrefoff := d.getr(onscan-1, phase=1, nfeed=1, nif=nif)

     ratio_out := (vsigoff.data.arr - vrefoff.data.arr) / (vrefon.data.arr - vrefoff.data.arr)
     diffs := (vsigoff.data.arr - vrefoff.data.arr)
     resids := (vsigon.data.arr - vsigoff.data.arr) - (vrefon.data.arr - vrefoff.data.arr)

     print onscan, median(ratio_out[1,50:975]), median(ratio_out[2,50:975]), median(resids[1,50:975]), median(resids[2,50:975]), median(diffs[1,50:975]), median(diffs[2,50:975]) 

}