; IDL procedure to compute the frequency of the 
; strongest features in a spectrum. Median filters
; the data and lists all features with amplitude
; greater than 4dB. Useful when identifying the
; test tone in a data set. Data must be collected 
; with the spectral processor in switched total power 
; mode using /GBTdata/setups/total.pwr.
;
pro Ttone, onscan, record

if record EQ 0 then begin
  print, "Record must be > 0. Setting record = 1..."
  record = 1
endif

!Quiet = 1  ; suppress compiler messages

 scanon = strtrim(onscan,2)

 for i = 1, 3, 1 do begin 
   fxbopen,t1,scanon,i,hdr
   extname = fxpar(hdr,'EXTNAME')
   name = strtrim(extname,2) 
     if !err lt 0 then extname = ''

   if name EQ 'DATA' then begin 
     last_record = fxpar(hdr, "NAXIS2")
     if record GT last_record then begin
       print, "Displaying data from record ", last_record
       record = last_record
     endif
     fxbread,t1,ondata,"DATA", record
     source = fxpar(hdr, 'OBJECT')
   endif

   if name EQ 'RECEIVER' then begin
     fxbread,t1,obsfreq,"OBSFREQ"
     fxbread,t1,rfside,"RFSIDE"
     fxbread,t1,ifside,"IFSIDE" 
     fxbread,t1,freqres,"FREQRES"
     fxbread,t1,bw,"BANDWD"
   endif
   fxbclose, t1
 endfor
 
 rawsz = size(ondata)
 nchans = rawsz(1)
 nphases = rawsz(2)
 nifs = rawsz(3)

 PRINT, "Channels = ", nchans
 PRINT, "Phases = ", nphases
 PRINT, "IFs = ", nifs

 cf = obsfreq(0) / 1e6
 PRINT, "Center Frequency = ", cf, " MHz"
 PRINT, "Bandwidth = ", bw(0) / 1e6, " MHz"


; rfsign = 1.0
; if (rfside(0) EQ 1) then rfsign = -rfsign
; if (ifside(0) EQ 1) then rfsign = -rfsign

 rfsign = -1.0
 f = dindgen(rawsz(1)) - (rawsz(1)/2)
 f = f * freqres(0) * rfsign + obsfreq(0)
 f = f / 1e6  ; units of MHz

; irow0 = ondata(*,0,0)
; imed0 = MEDIAN(irow0,13)
; irow0 = 10.0*ALOG10(irow0/imed0)

 irow1 = ondata(*,0,1)
 imed1 = MEDIAN(irow1,13)
 irow1 = 10.0*ALOG10(irow1/imed1)

; max0 = MAX(irow0)
 max1 = MAX(irow1)
; print, "Maximum values: ", max0, max1
 print, "Maximum value: ", max1


 FOR i = 0, nchans-1, 1 DO BEGIN
;   IF(irow0(i) GT 4.0) THEN PRINT, i, f(i), irow0(i)
   IF(irow1(i) GT 2.0) THEN PRINT, i, f(i), irow1(i)
 ENDFOR

 fmin = MIN(f)
 fmax = MAX(f)

; WINDOW, 0
; PLOT,  f, irow0, Xstyle=1, Ystyle=1, Xrange=[fmin,fmax], $
;        Yrange=[-1,2*max0], Xticks=4, Xminor=5, Yticks=10, Yminor=2, $
;        Ytitle='Normalized Power, dB', Xtitle='Frequency, MHz'

;SET_PLOT, 'PS'
 
 WINDOW, 1
 PLOT,  f, irow1, Xstyle=1, Ystyle=1, Xrange=[fmin,fmax], $
        Yrange=[-1,max1], Xticks=4, Xminor=5, Yticks=10, Yminor=2, $
        Ytitle='Normalized Power, dB', Xtitle='Frequency, MHz'

;device, /close
;file = 'idl.' + STRLOWCASE(!D.NAME)
;cmd = 'lpr -Plp ' + file
;spawn, cmd
;SET_PLOT, 'X

END