# This glish module retrieves autocorellation data from a autocorrelation
# spectrometer GBT FITS file and frequency information from the corresponding
# LO1A file and computes and displays the spectrum of the system temperature as
# a function of sky frequency.  It uses the FITS retrieval routine in the
# module 'get_fits', the autocorrelation correction functions from the
# 'VanVleck' module, and the FFT routine from the Numeric FFT module.  The FITS
# files are assumed to reside in the directory /home/gbtdata under the project
# name directory.  The AIPS++ plotting tool, 'pgplotter' which uses the glish
# interface to the PGPLOT library is used for plotting.  The AIPS++ FFT tool
# is also used.

# open FITS file access
include 'get_fits.g'
# open the plotting tool
if (!is_defined('pg')) {
    include 'pgplotter.g'
    pg := pgplotter();
    include 'fftserver.g'
    fft := fftserver();
}

# The arguments are the project name, e.g., 'AGBT02A_069', the FITS file
# name of the spectrometer data to be plotted, e.g.,
# '2003_05_24_20:44:16A.fits', and the calibration noise source values for each
# receiver channel, e.g., [1.7, 1.7].  The noise values are defaulted to 1.0.
# We assume that the only phase switching in the data is cal no, cal off.

func plot_ac_tsys ( project, file, tcal=[1.0, 1.0] )
{
    # assemble the directory paths and retrieve the spectral processor and
    # LO1A data records.  See the 'get_fits' module for the details of
    # the directory structure and a routine for listing a FITS file's
    # contents
    base_path := spaste('/home/gbtdata/', project, '/')
    acft := get_fits(spaste(base_path, 'Spectrometer/', file))
    # the spectrometer FITS file name has a bank designator suffix that needs
    # to be removed for use by other device file names
    base_file := get_base_file_name(file);
    loft := get_fits(spaste(base_path, 'LO1A/', base_file))

    # get the first mixer's sideband string
    rf_sideband := loft.main_header.SIDEBAND
    # compute the observed center frequency from the LO frequency at the
    # beginning of the scan, the first IF center frequency, and the sideband
    # and convert from Hz to MHz
    if (rf_sideband == 'LOWER') {
        ctr_freq := (loft.LO1TBL.Columns.LO1FREQ[1] -
                     loft.main_header.IFFREQ) * 1e-6
    } else {
        ctr_freq := (loft.LO1TBL.Columns.LO1FREQ[1] +
                     loft.main_header.IFFREQ) * 1e-6
    }
    # get the spectrum bandwidth in MHz.
    bandwidth := acft.PORT.Columns.BANDWDTH[1] * 1e-6
    # get the dimensions of the spectrum data array and give them physical
    # meaning
    data_shape := acft.DATA.Columns.DATA::shape
    num_records := data_shape[4] # number of integrations in the scan
    num_rcvrs := data_shape[3]   # number of IF channels into spectral processor
    num_phases := data_shape[2]  # switching phases, i.e., cal on, cal off
    num_chan := data_shape[1]    # number of frequency bins in each spectrum
    # compute the array of frequency offsets of the spectral channels with
    # respect to the reference channel, num_chan / 2, where the channel
    # numbers run from 1 to num_chan.  The number of frequency
    # channels after Fourier transforming the ACF is equal to the number of
    # ACF lags.
    ch_spacing := bandwidth / num_chan
    freq_offset := ([1:num_chan] - (num_chan / 2)) * ch_spacing
    # the GBT has three mixer stages before the spectral processor input,
    # but the second and third mixer conversions are both always lower sideband.
    # The effective sideband of the autocorrelation spectrometer depends on its
    # bandwidth, since there is an extra lower-sideband conversion before the
    # 50- and 12.5-MHz samplers, and the samplers use one of the aliases
    # rather than baseband.  The following rules apply:
    inversion_flag := 0
    bwkhz := as_integer(bandwidth * 1e3 + 0.5)  # kHz
    if ((bwkhz == 12500) || (bwkhz == 800000)) {
        inversion_flag := 1
    }
    if (((rf_sideband == 'LOWER') && (!inversion_flag)) ||
        ((rf_sideband == 'UPPER') && (inversion_flag))) {
        freq := ctr_freq - freq_offset
    } else {
        freq := ctr_freq + freq_offset
    }
    # if the frequency array runs from higher to lower, flip it and set a
    # flag for flipping the spectrum arrays as well
    reverse_flag := 0
    if (freq[1] > freq[len(freq)]) {
        freq := freq[len(freq):1]
        reverse_flag := 1
    }
    # set spectrum array phases for the two cal states (1 = on, 0 = off)
    if (acft.STATE.Columns.CAL[1] == 1) {
        on_state := 2
        off_state := 1
    } else {
        on_state := 1
        off_state := 2
    }
    # get the number of sampler levels used to use in thw correlation
    # quantization (Van Vleck) corrections.
    levels := acft.PORT.Columns.LEVEL[1]
    # assign the plot and axis titles
    xttl := 'Frequency in MHz'
    yttl := 'System Temperature in Kelvins'
    ttl := spaste('Autocorrelation Spectrometer Tsys Plot  ', file)
    # set the color of the first spectrum plot, the Y position of the
    # text legend, and the overplot flag
    color_name := ['White', 'Red  ', 'Green', 'Blue ', 'Aqua ']
    color := 2   # red
    txty := 1.5
    newplot := T
    pg.clear()
    # process the spectra for each receiver
    for (rcvr_num in [1:num_rcvrs]) {
        # if not enough cal values have been specified, default it to 1.0
        if (len(tcal) < rcvr_num) {
            tcal[rcvr_num] := 1.0
	}
        # create empty data arrays for the cal on and cal of spectrum
        # accumulations
        cal_on := array(0.0, num_chan)
        cal_off := array(0.0, num_chan)
        # accumulate the spectra from all integrations in the scan
        for (rec_num in [1:num_records]) {
            cal_on +:= get_spectrum(acft.DATA.Columns.DATA[ , rcvr_num,
						on_state, rec_num], levels)
            cal_off +:= get_spectrum(acft.DATA.Columns.DATA[ , rcvr_num,
						off_state, rec_num], levels)
	}
        # renormalize the spectra
        cal_on := cal_on / num_records
        cal_off := cal_off / num_records
	# do a Hanning convolution
        cal_on[2:(len(cal_on) - 1)] := 0.5 * (cal_on[2:(len(cal_on) - 1)] +
				        0.5 * (cal_on[1:(len(cal_on) - 2)] +
                                           cal_on[3:len(cal_on)]))
        cal_off[2:(len(cal_off) - 1)] := 0.5 * (cal_off[2:(len(cal_off) - 1)] +
				        0.5 * (cal_off[1:(len(cal_off) - 2)] +
                                           cal_off[3:len(cal_off)]))
        # reverse the spectra, if necessary
        if (reverse_flag) {
            cal_on := cal_on[len(tsys):1]
            cal_off := cal_off[len(tsys):1]
	}
        # compute the channel by channel system temperature and clip
        # out-of-range values
        tsys := tcal[rcvr_num] * cal_off / (cal_on - cal_off)
        # draw the spectrum
        pg.plotxy(freq, tsys, title=ttl, xtitle=xttl, ytitle=yttl,
		  linecolor=color, newplot=newplot)
	pg.setyscale(0.0, 30.0)
        # set the legend text X position
        txtx :=  freq[1] + 30 * ch_spacing
        # write the legend on the plot
        pg.text(txtx, txty, spaste(color_name[color], ': rcvr ', rcvr_num))
        # increment color and legend position
        color +:= 1       # 3 = green, 4 = blue, 5 = blue-green
        txty +:= 1.5
	newplot := F
    }
}

# This function removes the GBT spectrometer bank suffix from the FITS file
# name.

func get_base_file_name ( file )
{
    return spaste(split(file, 'ABCD'));
}

# open the quantization and linearization correction client
vv := client('/users/rfisher/Applications/GlishClients/VanVleck/vanvleck_client');

# This routine applies the quantization and linearization corrections to the
# autocorrelation function and returns the amplitude of the Fourier transform
# of the ACF

func get_spectrum ( acf, level )
{
    # 'optimum' is the maximum sensitivity threshold level for the ampler
    if (level == 3) {
	optimum := 0.5405
    } else {
	# the 9-level ACF received from the GBT spectrometer
	# is a factor of 16 smaller than the quantization mathematics assume
	acf := as_float(16.0 * acf)
	optimum := 3.401
    }
    # get the corrected ACF values and the threshold level and power
    # linerization correction derived from the zero lag value.  The returned
    # ACF is normalized to a zero-lag value of 1.
    cor_vals := vv->get_corrected_lags([acf=acf, levels=level]);
    # print some of this info
    zero_lag := acf[1]; zero_lag::print.precision := 4;
    frac := acf[1] / optimum; frac::print.precision := 4;
    thresh := cor_vals.threshold; thresh::print.precision := 4;
#    print spaste(level, '-level zero_lag:'), zero_lag, 'approx.',
#		 frac, 'of optimum, v_thresh:',thresh;
    # apply the linearization correction
    acf := cor_vals.acf * cor_vals.pwr_correction;
    # create a symmetric array out of the ACF by mirror image around its
    # zero-lag value
    dacf := array(0.0, 2 * length(acf));
    dacf[(length(acf) + 1):(2 * length(acf))] := acf;
    dacf[2:length(acf)] := acf[length(acf):2];
    dacf[1] := dacf[2];
    # do the Fourier transform which, because the input is real rather than
    # complex data, is redundant in its two halves.  Slice off only one half
    # to get an array running from zero+one channel to maximum frequency,
    # i.e., the DC channel (length(acf) + 1) is dropped.
    spec := fft.realtocomplexfft(dacf)[length(acf):1];
    # compute and return the complex vector amplitudes.
    return sqrt(real(spec)^2 + imag(spec)^2)
}

# To test, uncomment the plot call below and run :
# include 'plot_actsys.g'
print 'The quantization corrections for 40 spectra takes about 14 seconds:'
plot_ac_tsys('AGBT02A_069', '2003_05_24_20:44:16A.fits', [1.5, 1.5])