print, 'This program produces a time-frequency plot of DWE data files'
print, '                                   '
print, '                                   '
print, 'Output Medium:                     '
print, '                                   '
print, '               output.ps         (0)'
print, '               output.eps        (1)'
print, '               Tektronix-Screen  (2)'
print, '               X-Screen          (3)'
print, '               Windows           (4)'
print, '               terminate program (5)'
read, medium

case medium of
  0: begin
       set_plot, 'PS'
       device, filename='output.ps', $
;         bits = 8, /color, $
         /landscape,  $
           xoffset=2.5, yoffset=26., $
           xsize=24., ysize=16.
;         /portrait, $
;           xoffset=2., yoffset=3., $
;           xsize=16., ysize=24.
     end
  1: begin
       set_plot, 'PS'
       device, filename='output.eps', $
;         bits = 8, /color, $
         /encapsulated, $
         xsize=24., ysize=16.
     end
  2: begin
       set_plot, 'TEK'
     end
  3: begin
       set_plot, 'X'
;       window, 0, xsize=650.*210/297, ysize=650.
;       !p.color=255
;       !p.background=0
     end
  4: begin
       window, 0, xsize=600., ysize=400.
;       !p.color=255
;       !p.background=0
     end
  else: begin
          print, 'Program terminated'
          stop
        end
endcase

; --- Gray scale ---
;RED   = findgen(256)
;GREEN = findgen(256)
;BLUE  = findgen(256)
;tvlct, RED, GREEN, BLUE

dummy = ' '

print, 'Read data? (y/n)'
read, dummy
if dummy eq 'y' then begin
  readflag = 1
  openr, 1,'carriera.dat', error = err       ; t_start and t_end determination
  if err ne 0 then begin
    print, !err_string
    stop
  endif
  datc = dblarr(2)
  readf, 1, datc
  t_start = datc[0]
  while not eof(1) do readf, 1, datc
  close, 1
  t_end = datc[0]

  print, 'Time scale:'
  print, ' '
  print, '    t_0 = 0               (1)'
  print, '    Start of test = 0     (2)'
  print, '    scale with given time (3)'
  read, tscale

  case tscale of
    1: begin
         print, 'Enter moment of t_0 as recorded from HK3 (in hours)'
         read, t_sub
         xtit = '!18t - t!D0!N!X [min]'
       end
    2: begin
         t_sub = t_start
         xtit = 'TEST ELAPSED TIME [min] (T=0 <=> ' + string(floor(t_start), format = '(I0)') + ':' + string((t_start - floor(t_start)) * 60.0, format = '(D4.1)') + ' SCET)'
       end
    else: begin
            t_sub = 0.0
            xtit = 'SCET [min]'
          end
  endcase

  t_start = t_start - t_sub
  t_end = t_end - t_sub
endif else readflag = 0

print, 'Enter time resolution (in minutes)'
read, timeres
;timeres = 1.0

tname = ' '
tdate = ' '
if readflag eq 1 then begin
  print, 'Enter test name'
  read, tname
  print, 'Enter date of test'
  read, tdate
endif

print, 'Enter minimum frequency-range (-1 for range 0-4 Hz)'
read, minfreq
;minfreq = 0.0
;maxfreq = 4.0
;maxfreq = 1.0
if minfreq ne -1 then begin
  print, 'Enter maximum frequency-range'
  read, maxfreq
;  minfreq = 0.001
;  maxfreq = 1.0
endif else begin
  minfreq = 0.0
  maxfreq = 4.0
endelse

ylogflag = 0
;print, 'Use logarithmic frequency axis? (y/n)'
;read, dummy
;dummy = 'n'
;if dummy eq 'y' then ylogflag = 1

print, 'Enter minimum time-range (-1 for full range)'
read, mintime
;mintime = -1
if mintime ne -1 then begin
   print, 'Enter maximum time-range'
   read, maxtime
endif else begin
   mintime = t_start * 60.0
   maxtime = t_end * 60.0
endelse

if readflag eq 1 then begin
  smflag = 0
  print, 'Use residuals (orig. data - 10-sec-smooth)? (y/n)'
  read, dummy
  if dummy eq 'y' then smflag = 1

;  print, 'File:'
;  print, ' '
;  print, '        NCO-A         (1)'
;  print, '        NCO-B         (2)'
;  print, '        AGC-A         (3)'
;  print, '        AGC-B         (4)'
;  read, fname_no
  fname_no = 1
  case fname_no of
    1: begin
         file = 'NCO-A'
         fname = 'ncoa.dat'
         ytit = 'f!DR!N (A) [Hz]'
       end
    2: begin
         file = 'NCO-B'
         fname = 'ncob.dat'
         ytit = 'f!DR!N (B) [Hz]'
       end
    3: begin
         file = 'AGC-A'
         fname = 'agca.dat'
         ytit = 'AGC-A [dBm]'
       end
    4: begin
         file = 'AGC-B'
         fname = 'agcb.dat'
         ytit = 'AGC-B [Hz]'
       end
    else: begin
            print, 'Not a valid number, program terminated'
            stop
          end
  endcase

;  tit = 'Spectrum: ' + tname + ', ' + tdate + ', ' + file + $
;    ', Time Resolution: ' + string(timeres, format='(i2)') + ' min'
  if smflag eq 1 then $
    ytit = ytit + ' (Residuals)'
endif
tit = 'Spectrum: ' + tname + ', ' + tdate + ', ' + file + $
  ', Time Resolution: ' + string(timeres, format='(d4.1)') + ' min'

xyouts, 0.5, 0.5, /normal, '!5'  ; setzt den Font

!x.thick = 3.0
!y.thick = 3.0
!p.charsize = 1.5
!p.charthick = 2.0

datplotflag = 0
;print, 'Plot data as well? (y/n)'
;read, dummy
dummy = 'y'
if dummy eq 'y' then begin
  datplotflag = 1
  !p.multi = [0, 1, 3, 0, 0]

;  print, 'Enter minimum Y-range for data-plot (-1 for auto detection)'
;  read, ymin
  ymin = -20.
;  if ymin ne -1 then begin
;     print, 'Enter maximum Y-range'
;     read, ymax
;  endif else begin
;     ymin = 0.0
;     ymax = 0.0
;  endelse
  ymax = 20.
  print, 'Connect data points by lines? (y/n)'
  read, dummy
;  dummy = 'n'
  if dummy eq 'n' then psy = 3 else psy = 0
endif else !p.multi = 0

;t_index0 = 0.0
;read, t_index0, prompt = 'Enter time offset of first data point for FFT (in minutes): '
;index = long(floor(t_index0 * 480.0))
index = 0

if readflag eq 1 then begin
  nolines = long(0.0)               ; count lines
  openr, 1, fname
  while not eof(1) do begin
     readf, 1, dummy
     nolines = nolines + 1
  endwhile
  close, 1

  dat = dblarr(2,nolines)           ; read data
  openr, 1, fname
  readf, 1, dat
  close, 1
  dat[0,*] = (dat[0,*] - t_sub) * 60.0  ; <=> minutes
  if (smflag eq 1) then dat[1,*] = dat[1,*] - smooth(dat[1,*], 81)
endif

nop = floor(timeres * 480.0)  ; no. of points = res. * 60 sec/min * 8 samples/sec
; nob = ceil(double(nolines) / double(nop))   ; no. of time bins + 1 for the last bin
nob = floor(double(nolines-index) / double(nop))   ; no. of time bins

ampl = dcomplexarr(nob, nop)
time = dblarr(nob)
freq = dindgen(nop) * 8.0 / (nop-1)

for i=long(0.0), nob-1 do begin

  ampl[i,*] = fft(dat[1, index:index+nop-1], 1, /double)     ; array !!
  time[i] = (dat[0,index] + dat[0,index+nop-1]) / 2.0

  index = index + nop
  if index + nop ge nolines then i = nob-1

endfor

powdens = 2.0 * abs(ampl) * abs(ampl) / nop / 8.0
pdeq0 = where(powdens le 0.0, count)
if count ne 0 then powdens(pdeq0) = 1.e-20

if fname_no le 2 then begin
  lev = [1., 2., 5., 1.e1, 2.e1, 5.e1, 1.e2, 2.e2, 5.e2, 1.e3, 2.e3, 5.e3, 1.e4, 3.e4, 5.e4, 1.e5, 2.e5, 5.e5, 1.e6]
  noiselev = 10.0
endif else begin
  lev = [1.e-2, 2.e-2, 5.e-2, 1.e-1, 2.e-1, 5.e-1, 1.e0, 2.e0, 5.e0, 1.e1, 2.e1, 5.e1, 1.e2, 2.e2, 5.e2, 1.e3, 2.e3, 5.e3]
  noiselev = 0.1
endelse

lev = alog10(lev)
noiselev = alog10(noiselev)

if datplotflag eq 1 then begin
  !p.region   = [0.0, 0.65, 0.9, 0.95]
  !p.position = [0.1, 0.65, 0.9, 0.95]
  dummy = ' '
  !x.tickname = replicate(' ', 30)
endif else begin
  dummy = xtit
  !x.tickname = replicate('', 30)
endelse

;contour, powdens, time, freq, $
contour, alog10(powdens), time, freq, $
;contour, min_curve_surf(powdens, xvalues = time, yvalues = freq), time, freq, $
;  /closed, $   ; sieht scheisse aus
;  /fill, $    ; nur bei Farbe
;  c_color = 200-200*indgen(n_elements(lev))/(n_elements(lev)-1), $
  title = tit, $
;  c_labels = replicate(0,30), $
  ylog = ylogflag, $
  xrange = [mintime, maxtime], $
  xstyle = 1, $
  xtitle = dummy, $
  yrange = [minfreq, maxfreq], $
;  ystyle = 1, $
  ystyle = 9, $
  ytitle = 'f!DFFT!N [Hz]', $
  levels = lev, $
  c_linestyle = (lev lt noiselev)

axis, /yaxis, ystyle = 1, $
  yrange = [minfreq/0.3663194, maxfreq/0.3663194], $
  ytitle = 'Multiples of f!Dosc!N'

if datplotflag eq 1 then begin
  !p.region   = [0.0, 0.35, 0.9, 0.65]
  !p.position = [0.1, 0.35, 0.9, 0.65]

  if fname_no le 2 then ysty = 9 else ysty = 1
  plot, dat[0,*], dat[1,*], $
    psym = psy, $
    xrange = [mintime, maxtime], $
    xstyle = 1, $
;    xtitle = xtit, $
    xtitle = ' ', $
    yrange = [ymin, ymax], $
    ystyle = ysty, $
;    ytitle = 'f!DR!N [Hz]', $
    ytitle = ytit, $
;    xtickname = replicate('', 30)
    xtickname = replicate(' ', 30)

  if fname_no le 2 then begin
    axis, /yaxis, ystyle = 1, $
      yrange = [ymin / 6.8, ymax / 6.8], $
      ytitle = 'v!Dr!N [m/s]'
  endif
endif

; additional AGC-Plot for PRT evaluation
;nolines = long(0.0)               ; count lines
;openr, 1, 'agca.dat'
;while not eof(1) do begin
;  readf, 1, dummy
;  nolines = nolines + 1
;endwhile
;close, 1

datagc = dblarr(2,nolines)           ; read data
openr, 1, 'agca.dat'
readf, 1, datagc
close, 1
datagc[0,*] = (datagc[0,*] - t_sub) * 60.0  ; <=> minutes

!p.region   = [0.0, 0.05, 0.9, 0.35]
!p.position = [0.1, 0.05, 0.9, 0.35]
plot, datagc[0,*], datagc[1,*], $
  psym = psy, $
  xrange = [mintime, maxtime], $
  xstyle = 1, $
  xtitle = xtit, $
  yrange = [0,0], $
  ystyle = 1, $
  ytitle = 'AGC [dBm]', $
  xtickname = replicate('', 30)

;; gestrichelte Linien fuer Kalibration
;cal = 549
;oplot, [!x.crange[0], !x.crange[1]], [cal, cal], $
;  linestyle = 2, thick = 0.5
;xyouts, !x.crange[1], cal, '-130 dBm', charsize = .6
;cal = 351
;oplot, [!x.crange[0], !x.crange[1]], [cal, cal], $
;  linestyle = 2, thick = 0.5
;xyouts, !x.crange[1], cal, '-120 dBm', charsize = .6
;cal = 121
;oplot, [!x.crange[0], !x.crange[1]], [cal, cal], $
;  linestyle = 2, thick = 0.5
;xyouts, !x.crange[1], cal, '-108 dBm', charsize = .6

if medium eq 0 or medium eq 1 then device, /close

finished:
!x.range = 0
!x.style = 0
!x.tickname = 0
!y.style = 0
!x.thick = 1.0
!y.range = 0
!y.thick = 1.0
!p.multi = 0
!p.psym = 0
!p.charsize = 1.0
!p.charthick = 1.0
!p.title = ''
!p.region = 0
!p.position = 0

end