; This IDL 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.  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_actsys, 'AGBT02A_069', '2003_05_24_20:44:16A.fits', [1.5, 1.5]

; This code does not include quantization or lineariztion corrections to the
; autocorrelation function

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

function get_base_file_name, file
    tpos = 100
    if (STRPOS(file, 'A') GT 0) then begin
	tpos = STRPOS(file, 'A')
    endif
    if (STRPOS(file, 'B') GT 0) then begin
	tpos = STRPOS(file, 'B')
    endif
    if (STRPOS(file, 'D') GT 0) then begin
	tpos = STRPOS(file, 'C')
    endif
    if (STRPOS(file, 'D') GT 0) then begin
	tpos = STRPOS(file, 'D')
    endif
    if tpos LT 100 then begin
	pfx = STRMID(file, 0, tpos)
	sfx = STRMID(file, tpos + 1, 100)
	return, pfx + sfx
    endif else begin
	return, file
    endelse
end

function power_correction, zeroLag, lev

    if (lev EQ 3) then begin
        erfData = (1.0 - zeroLag)
        ipl = erfData * erfData - 0.5625
        plsq = ipl * ipl
        ipl = erfData * (1.591863138 + (ipl * (-2.442326820)) + $
                                      (plsq * 0.37153461)) / (1.467751692 + $
                                               (ipl * (-3.013136362)) + plsq)
        return, 0.3745443672 / (2.0 * ipl * ipl)
    endif else begin
        zeroLag = zeroLag / 16.0;
        coef0 =  -0.03241744594;
        coef1 =   4.939640303;
        coef2 =  -5.751574913;
        coef3 =  34.83143031;
        coef4 = -78.66637472;
        coef5 = 213.7108496;
        coef6 =-317.1011469;
        coef7 = 245.8618017;

        return, ((((((coef7 * zeroLag + coef6) * zeroLag + $
                    coef5) * zeroLag + coef4) * zeroLag + $
                  coef3) * zeroLag + coef2) * zeroLag + $
                coef1) * zeroLag + coef0
    endelse
end

; This function does a Fourier transform on the autocorrelation function and
; returns the spectrum.  This is also where the quantization (Van Vleck)
; and lineraization corrections sould be applied, when implemented.

function get_spectrum, acf, levels

    ; 'optimum' is the maximum sensitivity threshold level for the ampler
    if (levels EQ 3) then begin
	optimum = 0.5405
        ; 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.
	pwr_cor = power_correction(acf(0), 3)
; do the 3-level quantization correction here but, for now, just normalize
; the ACF as would be done in the correction
	acf = acf / acf(0)
    endif else begin
	; the 9-level ACF received from the GBT spectrometer
	; is a factor of 16 smaller than the quantization mathematics assume
	acf = 16.0 * acf
	pwr_cor = power_correction(acf(0), 9)
	optimum = 3.401
; do the 3-level quantization correction here but, for now, just normalize
; the ACF as would be done in the correction
	acf = acf / acf(0)
    endelse
    ; apply the linearization correction
    acf = pwr_cor * acf
    ; create a symmetric array out of the ACF by mirror image around its
    ; zero-lag value
    sizes = SIZE(acf)
    len_acf = sizes(1)
    dacf = FLTARR(2 * len_acf)
    dacf(len_acf:*) = acf
    dacf(1:(len_acf - 1)) = REVERSE(acf(1:*))
    dacf(0) = dacf(1);
    ; do the Fourier transform
    spec = FFT(dacf)
    ; slice off a non-redundant half
    spec = spec(1:len_acf)
    ; return the amplitude of the complex vectors
    return, SQRT(spec * CONJ(spec))
end

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

pro plot_actsys, project, file, tcal

    ; 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 = '/home/gbtdata/' + project + '/'
    acfn = 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)
    lofn = base_path + 'LO1A/' + base_file
    ; check for the existence of the FITS files
    GET_LUN, fd
    OPENR, fd, acfn
    stat = FSTAT(fd)
    if (stat.OPEN EQ 0) then begin
	print, 'Cannot open FITS file: ', acfn
	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 PORT HDU
    fxbopen, t1, acfn, 1, ehdr
    ; get the spectrum bandwidth in MHz.
    fxbread, t1, bw, 'BANDWDTH'
    bandwidth = bw(0)
    ; get the number of sampler levels used to use in the correlation
    ; quantization (Van Vleck) corrections.
    fxbread, t1, level, 'LEVEL'
    levels = level(0)
    fxbclose, t1

    ; open the STATE HDU
    fxbopen, t1, acfn, 2, ehdr
    ; get the cal states
    fxbread, t1, cal_state, 'CAL'
    bandwidth = bw(0) * 1e-6
    fxbclose, t1

    ; open the DATA HDU
    fxbopen, t1, acfn, 5, ehdr
    ; get the autocorrelation array
    fxbread, t1, acf, 'DATA'
    fxbclose, t1

    ; get the dimensions of the autocorrelation data array and give them
    ; physical meaning
    sizes = SIZE(acf)
    num_chan = sizes(1)		; number of frequency bins in each spectrum
    num_rcvrs = sizes(2)	; number of IF channels into spectral processor
    num_phases = sizes(3)	; switching phases, i.e., cal on, cal off
    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 - 1, where the channel
    ; numbers run from 0 to num_chan - 1.  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 = (INDGEN(num_chan) - (num_chan / 2) + 1) * 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 = ROUND(bandwidth * 1e3)  ; kHz
    if (bwkhz EQ 12500) OR (bwkhz EQ 800000) then inversion_flag = 1
    if ((rf_sideband EQ 'LOWER') AND (inversion_flag EQ 0)) OR $
       ((rf_sideband EQ 'UPPER') AND (inversion_flag 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)
    if cal_state(0) EQ 1 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 + $
		      get_spectrum(acf(*, rcvr_num, on_state, rec_num), levels)
            cal_off = cal_off + $
		      get_spectrum(acf(*, rcvr_num, off_state, rec_num), levels)
	endfor
        ; renormalize the spectra
        cal_on = cal_on / num_records
        cal_off = cal_off / num_records
	; do a Hanning convolution
        cal_on[1:(num_chan - 2)] = 0.5 * (cal_on[1:(num_chan - 2)] + 0.5 * $
					  (cal_on[0:(num_chan - 3)] + $
                                           cal_on[2:(num_chan - 1)]))
        cal_off[1:(num_chan - 2)] = 0.5 * (cal_off[1:(num_chan - 2)] + 0.5 * $
					  (cal_off[0:(num_chan - 3)] + $
                                           cal_off[2:(num_chan - 1)]))
        ; reverse the spectra, if necessary
        if reverse_flag EQ 1 then begin
	    cal_on = REVERSE(cal_on)
	    cal_off = REVERSE(cal_off)
	endif
        ; 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
	if rcvr_num EQ 0 then begin
;	    PLOT, freq, cal_off
;	    OPLOT, freq, cal_on
	    !X.STYLE = 1
	    PLOT, freq, tsys, YRANGE=[0, 30], $
		XRANGE=[freq(0), freq(num_chan - 1)], $
		TITLE='Autocorrelation Spectrometer Tsys Plot  ' + file, $
		XTITLE='Frequency in MHz', $
		YTITLE='System Temperature in Kelvins'
	endif else begin
	    OPLOT, freq, tsys, COLOR=255
	    dummy = 0
	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
        ; increment color and legend position
        txty = txty + 1.5
    endfor
end