# This glish script retrieves spectrum data from a spectral processor 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 script
# 'get_fits'.  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.

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

# The arguments are the project name, e.g., 'AGBT02A_069_01', the FITS file
# name of the data to be plotted, e.g., '2003_06_17_09:22:34.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_sp_tsys ( project, file, tcal=[1.0, 1.0] )
{
    # assemble the directory paths and retrieve the spectral processor and
    # LO1A data dictionaries.  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, '/')
    spft := get_fits(spaste(base_path, 'SpectralProcessor/', file))
    loft := get_fits(spaste(base_path, 'LO1A/', 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 spacing of the spectral processor spectrum in MHz
    ch_spacing := spft.RECEIVER.Columns.FREQRES[1] * 1e-6
    # get the spectral processor input sideband: 0 = upper so that increasing
    # output spectrum channel number corresponds to frequency at the spectral
    # processor input.  if_sideband = 1 is the reverse.
    if_sideband := spft.RECEIVER.Columns.IFSIDE[1]
    # get the dimensions of the spectrum data array and give them physical
    # meaning
    data_shape := spft.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 frquency offsets of the spectral channels with
    # respect to the reference channel, num_chan / 2, where the channel
    # numbers run from 0 to num_chan - 1.
    freq_offset := ([0:(num_chan - 1)] - (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
    # so the following frequency offset to RF frequency rules apply:
    if (((rf_sideband == 'LOWER') && (if_sideband == 0)) ||
        ((rf_sideband == 'UPPER') && (if_sideband == 1))) {
        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 (spft.STATE.Columns.CAL[1] == 0) {
        on_state := 2
        off_state := 1
    } else {
        on_state := 1
        off_state := 2
    }
    # assign the plot and axis titles
    xttl := 'Frequency in MHz'
    yttl := 'System Temperature in Kelvins'
    ttl := spaste('Spectral Processor 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[rvcr_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 +:= spft.DATA.Columns.DATA[ , on_state, rcvr_num, rec_num]
            cal_off +:= spft.DATA.Columns.DATA[ , off_state, rcvr_num, rec_num]
	}
        # renormalize the spectra
        cal_on := cal_on / num_records
        cal_off := cal_off / num_records
        # compute the channel by channel system temperature
        tsys := tcal[rcvr_num] * cal_off / (cal_on - cal_off)
        # reverse the spectrum, if necessary
        if (reverse_flag) {
            tsys := tsys[len(tsys):1]
	}
        # 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
    }
}

# To test uncomment the plot call below and run :
# import plot_tsys
#   or
# reload(plot_tsys)
# plot_sp_tsys('AGBT02A_069_01', '2003_06_17_09:22:34.fits', [1.7, 1.7])