print, 'Ausgabemedium:                     '
print, '                                   '
print, '               output.ps        (0)'
print, '               Tektronix-Schirm (1)'
print, '               X-Schirm         (2)'
print, '               Windows          (3)'
print, '               Programmende     (4)'
read, medium

case medium of
  0: begin
       set_plot, 'PS'
       device, filename='output.ps', $
;         /landscape,  $
;           xoffset=2.5, yoffset=26., $
;           xsize=24., ysize=16.

           /portrait, $
            xoffset=2., yoffset=3., $
            xsize=16., ysize=24.
     end
  1: begin
       set_plot, 'TEK'
     end
  2: begin
       set_plot, 'X'
;       window, 0, xsize=650.*210/297, ysize=650.
;       !p.color=255
;       !p.background=0
     end
  3: begin
       window, 0, xsize=600., ysize=400.
;       !p.color=255
;       !p.background=0
     end
  else: begin
          print, 'Programmende'
          stop
        end
endcase

!x.range = 0

dummy = ' '

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

!p.multi = [0, 1, 2, 0, 0]      ; legt fest, wieviele plots auf 1 seite kommen

!p.psym = 3                      ; keine linien sondern punkte

;!p.charsize = 2.0
!p.charthick = 2.0

!x.thick = 3.0
!y.thick = 3.0

;print, 'Read data ?? (y/n)'
;read, dummy
;if dummy eq 'y' then begin
  openr, 1,'ncoa.dat'
  ln = long(0.0)
  while not eof(1) do begin
     readf, 1, dummy
     ln = ln + 1
  endwhile
  close, 1

  dat = dblarr(2,ln)
  dat_sm = dblarr(2,ln)
  openr, 1,'ncoa.dat'
  readf, 1, dat
  close, 1
  t_start = dat(0,0)
;  dat(0,*) = (dat(0,*) - t_start) * 60.0  ; <=> Minuten
  dat(0,*) = (dat(0,*) - t_start) * 60.0
;endif

lb = 130.0
rb = 160.0

; 1. NCO-Plot, probe level
;np = 81
np = 9

dat_sm(0,*) = dat(0,*)
dat_sm(1,*) = smooth(dat(1,*), np)

datmean = 0.0
datanzahl = 0.0
for i=long(0.0), ln-1 do begin
   if ((dat(0,i) ge !x.crange(0)) and (dat(0,i) le !x.crange(1))) then begin
      datmean = datmean + dat(1,i)
      datanzahl = datanzahl + 1.0
   endif
endfor
datmean = datmean / datanzahl
rms = 0.0
for i=long(0.0), ln-1 do begin
   if ((dat(0,i) ge !x.crange(0)) and (dat(0,i) le !x.crange(1))) then $
      rms = rms + (dat(1,i) - datmean) * (dat(1,i) - datmean)
endfor
rms = sqrt(rms / datanzahl)

datmean_sm = 0.0
datanzahl_sm = 0.0
for i=long(0.0), ln-1 do begin
   if ((dat_sm(0,i) ge !x.crange(0)) and (dat_sm(0,i) le !x.crange(1))) then begin
      datmean_sm = datmean_sm + dat_sm(1,i)
      datanzahl_sm = datanzahl_sm + 1.0
   endif
endfor
datmean_sm = datmean_sm / datanzahl_sm
rms_sm = 0.0
for i=long(0.0), ln-1 do begin
   if ((dat(0,i) ge !x.crange(0)) and (dat(0,i) le !x.crange(1))) then $
      rms_sm = rms_sm + (dat_sm(1,i) - datmean_sm) * (dat_sm(1,i) - datmean_sm)
endfor
rms_sm = sqrt(rms_sm / datanzahl_sm)

!p.region = [0.0, 0.5, 1.0, 1.0]
plot, dat(0,*), dat(1,*) - datmean_sm, $
;   position = [0.1, 0./18., 1.0, 1./18.], $
;   xtitle = 'TIME AFTER START OF TEST [min]', $
   xtitle = ' ', $
   xstyle = 1, $
   ytitle = 'f!DR!N [Hz]', $
   ystyle = 1, $
   yrange = [-10., 10.], $
   xrange = [lb, rb], $
   title = 'FM Probe Level Test: IST#1_CO#2 (13.8.96)'

;oplot, dat_sm(0,*), dat_sm(1,*) - datmean, psym = 0
oplot, dat(0,*), dat(1,*) - datmean, $
  nsum = 9, psym = 0

;xyouts, 132.0, 8.5, /data, $
;   '<f!DR!N> = ' + string(datmean, format='(d6.2)') + ' Hz', $
;   charsize = 1., charthick = 2.0
;xyouts, 132.0, 7.0, /data, $
;   'RMS = ' + string(rms, format='(d5.2)') + ' Hz', $
;   charsize = 1., charthick = 2.0
xyouts, 143.0, 8.5, /data, $
   'data smoothed over ' + string(np*0.125, format='(d6.3)') + ' sec', $
;   'data smoothed over 10 sec', $
   charsize = 1., charthick = 2.0
xyouts, 143.0, 7.0, /data, $
   'RMS (smoothed) = ' + string(rms_sm, format='(d5.2)') + ' Hz', $
   charsize = 1., charthick = 2.0

; 2. NCO-Plot, unit level

;if dummy eq 'y' then begin
  openr, 1, 'freq-ul_1s.dat'
  ln = long(0.0)
  while not eof(1) do begin
     readf, 1, dummy
     ln = ln + 1
  endwhile
  close, 1

  dat = dblarr(2,ln)
;  dat_sm = dblarr(2,ln)
  openr, 1,'freq-ul_1s.dat'
  readf, 1, dat
  close, 1
  t_start = dat(0,0)
;  dat(0,*) = (dat(0,*) - t_start) / 60.0  ; <=> Minuten
  dat(0,*) = (dat(0,*) - t_start) / 60.0
  dat(1,*) = dat(1,*) * 204.0
;endif

datmean = 0.0
datanzahl = 0.0
for i=long(0.0), ln-1 do begin
   if ((dat(0,i) ge !x.crange(0)) and (dat(0,i) le !x.crange(1))) then begin
      datmean = datmean + dat(1,i)
      datanzahl = datanzahl + 1.0
   endif
endfor
datmean = datmean / datanzahl
rms = 0.0
for i=long(0.0), ln-1 do begin
   if ((dat(0,i) ge !x.crange(0)) and (dat(0,i) le !x.crange(1))) then $
      rms = rms + (dat(1,i) - datmean) * (dat(1,i) - datmean)
endfor
rms = sqrt(rms / datanzahl)

!p.region = [0.0, 0.0, 1.0, 0.5]
plot, dat(0,*), dat(1,*) - datmean, $
;   position = [0.1, 0./18., 1.0, 1./18.], $
   xtitle = 'TIME AFTER START OF TEST [min]', $
   xstyle = 1, $
   ytitle = 'f!DR!N [Hz]', $
   ystyle = 1, $
   yrange = [-10.0, 10.0], $
   xrange = [lb, rb], $
   psym = 0, $
   title = 'FM Unit Level Test (6.3.96)'

;xyouts, 132.0, 8.5, /data, $
;   '<f!DR!N> = ' + string(datmean, format='(d6.2)') + ' Hz', $
;   charsize = 1., charthick = 2.0
;xyouts, 132.0, 7.0, /data, $
;   'RMS = ' + string(rms, format='(d5.2)') + ' Hz', $
;   charsize = 1., charthick = 2.0
xyouts, 143.0, 8.5, /data, $
;   'data smoothed over ' + string(np*0.125, format='(d6.3)') + ' sec', $
   'data smoothed over 1 sec', $
   charsize = 1., charthick = 2.0
xyouts, 143.0, 7.0, /data, $
   'RMS (smoothed) = ' + string(rms, format='(d5.2)') + ' Hz', $
   charsize = 1., charthick = 2.0

if medium eq 0 then device, /close

!x.range = 0
!x.thick = 1.0
!y.range = 0
!y.style= 1
!y.thick = 1.0
!p.multi = 0
!p.psym = 0
!p.charsize = 1.0
!p.charthick = 1.0
!p.title = ''
!p.region = 0

end