; This IDL module 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.  The FITS files are assumed to reside in the directory
; /home/gbtdata under the project name directory.

; To run from the IDL prompt:
; IDL> plot_tsys, 'AGBT02A_069_01', '2003_06_17_09:22:34.fits', [1.5, 1.5]

; 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].
; We assume that the only phase switching in the data is cal on, cal off.

pro plot_tsys, project, file, tcal

    ; assemble the directory paths
    base_path = '/home/gbtdata/' + project + '/'
    spfn = base_path + 'SpectralProcessor/' + file
    lofn = base_path + 'LO1A/' + file
    ; check for the existence of the FITS files
    GET_LUN, fd
    OPENR, fd, spfn
    stat = FSTAT(fd)
    if (stat.OPEN EQ 0) then begin
	print, 'Cannot open FITS file: ', spfn
	return
    endif
    CLOSE, fd
    GET_LUN, fd
    OPENR, fd, lofn
    stat = FSTAT(fd)
    if (stat.OPEN EQ 0) then begin
	print, 'Cannot open FITS file: ', lofn
	return
    endif
    CLOSE, fd

    ; get the LO FITS file main header
    mhdr = HEADFITS(lofn, EXTEN=0)
    ; get the first mixer's sideband string, and the first IF center frequency
    rf_sideband = fxpar(mhdr, 'SIDEBAND')
    if_freq = fxpar(mhdr, 'IFFREQ')
    ; get the LO frequency throughout the scan
    fxbopen, t1, lofn, 3, ehdr
    fxbread, t1, lo_freq, 'LO1FREQ'
    fxbclose, t1
    ; 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
    rf_sideband = STRTRIM(rf_sideband, 2)
    if (rf_sideband EQ 'LOWER') then begin
        ctr_freq = (lo_freq(0) - if_freq) * 1e-6
    endif else begin
        ctr_freq = (lo_freq(0) + if_freq) * 1e-6
    endelse
    ; open the RECEIVER HDU
    fxbopen, t1, spfn, 1, ehdr
    ; get the spacing of the spectral processor spectrum in MHz
    fxbread, t1, freq_res, 'FREQRES'
    ch_spacing = freq_res(0) * 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.
    fxbread, t1, if_side, 'IFSIDE'
    if_sideband = FIX(if_side(0))
    fxbclose, t1
    ; get the spectrum data array and give its dimensions physical meaning
    fxbopen, t1, spfn, 3, ehdr
    fxbread, t1, spec, 'DATA'
    fxbclose, t1
    sizes = SIZE(spec)
    num_chan = sizes(1)		; number of frequency bins in each spectrum
    num_phases = sizes(2)	; switching phases, i.e., cal on, cal off
    num_rcvrs = sizes(3)	; number of IF channels into spectral processor
    num_records = sizes(4)	; number of integrations in the scan
    ; 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 = (INDGEN(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
    ; so the following frequency offset to RF frequency rules apply:
    if (((rf_sideband EQ 'LOWER') AND (if_sideband EQ 0)) OR $
        ((rf_sideband EQ 'UPPER') AND (if_sideband EQ 1))) then begin
        freq = ctr_freq - freq_offset
    endif else begin
        freq = ctr_freq + freq_offset
    endelse
    ; 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(0) GT freq(num_chan - 1)) then begin
        freq = REVERSE(freq)
        reverse_flag = 1
    endif
    ; set spectrum array phases for the two cal states (1 = on, 0 = off)
    fxbopen, t1, spfn, 2, ehdr
    fxbread, t1, cal, 'CAL'
    fxbclose, t1
    if (cal(0) EQ 0) then begin
        on_state = 1
        off_state = 0
    endif else begin
        on_state = 0
        off_state = 1
    endelse
    ; set the Y position of the text legend
    txty = 1.5
    ; process the spectra for each receiver
    for rcvr_num = 0, (num_rcvrs - 1), 1 do begin
        ; create empty data arrays for the cal on and cal of spectrum
        ; accumulations
        cal_on = FLTARR(num_chan)
        cal_off = FLTARR(num_chan)
        ; accumulate the spectra from all integrations in the scan
        for rec_num = 0, (num_records - 1), 1 do begin
            cal_on = cal_on + spec(*, on_state, rcvr_num, rec_num)
            cal_off = cal_off + spec(*, off_state, rcvr_num, rec_num)
	endfor
        ; renormalize the spectra
        cal_on = cal_on / num_records
        cal_off = cal_off / num_records
        ; compute the channel by channel system temperature and clip
        ; out-of-range values
        tsys = tcal(rcvr_num) * cal_off / (cal_on - cal_off)
        ; reverse the spectrum, if necessary
        if reverse_flag EQ 1 then tsys = REVERSE(tsys)
        ; draw the spectrum
	if rcvr_num EQ 0 then begin
	    !X.STYLE = 1
	    PLOT, freq, tsys, YRANGE=[0, 30], $
		XRANGE=[freq(0), freq(num_chan - 1)], $
		TITLE='Spectral Processor Tsys Plot  ' + file, $
		XTITLE='Frequency in MHz', $
		YTITLE='System Temperature in Kelvins'
	endif else begin
	    OPLOT, freq, tsys, COLOR=255
	endelse
        ; set the legend text X position
        txtx =  freq[0] + 30 * ch_spacing
        ; write the legend on the plot
	if (rcvr_num EQ 0) then begin
	    clr = ' white'
	endif else begin
	    clr = ' red'
	endelse
	XYOUTS, txtx, txty, 'Rcvr ' + STRTRIM(STRING(rcvr_num + 1), 2) + clr
        txty = txty + 1.5
    endfor
end