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
medium = 0

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, 'PS'
	device, filename='output.eps', $
	/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
       set_plot, 'WIN'
;       window, 0, xsize=600., ysize=400.
;       !p.color=255
;       !p.background=0
     end
  else: begin
          print, 'Program terminated'
          stop
        end
endcase

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

xyouts, 0.5, 0.5, /normal, '!5' ; set font to font-no.5
 ; set font to font-no.5

dummy = ' '

;fn = 'recov.out'
fn = ' '
read, fn, prompt = 'Enter file name for ZWR: '

fflag = ' '
read, fflag, prompt = 'Enter zonal wind direction (p/r/0): '
if fflag eq 'p' then begin
  zwd = 1.
  yr = [-10., 140.]
endif
if fflag eq 'r' then begin
  zwd = -1.
  yr = [-140., 10.]
endif
if fflag eq '0' then begin
  zwd = 0.
  yr = [-10., 10.]
endif

abs_yrres = 0.
read, abs_yrres, prompt = 'Enter range for residual plot (e.g. 3): '
yrres = [-abs_yrres, abs_yrres]

ln = long(0.)
col = 16

openr, 1, fn
while not eof(1) do begin
  readf, 1, dummy
  ln = ln + 1
endwhile
close, 1

dat = dblarr(col, ln)

openr, 1, fn
readf, 1, dat
close, 1

!p.multi = [0,1,2,0,0]

!p.region = [0., 0.5, 1.0, 1.0]
plot, dat[0,*], dat[15,*], $
  position = [0.1, 0.5, 1.0, 1.0], $
;  xtitle = '!18t - t!D0!N!X [min]', $
  xtitle = ' ', $
  ytitle = 'Zonal Wind Speed [m/s]', $
;  yrange = [-10., 140.], $
  yrange = yr, $
  ystyle = 1, $
;  xrange = [0., 180.], $
  xstyle = 9, $
  nsum = 10, $
  xtickname = replicate(' ', 30)

oplot, dat[0,*], dat[15,*] + dat[14,*], $
  linestyle = 1

oplot, dat[0,*], dat[15,*] - dat[14,*], $
  linestyle = 1

for alt = 160, 0, -10 do begin
  if (alt le dat[1,0]) and (alt ge dat[1,ln-1]) then begin
    j = long(0.)
    while (not ((alt le dat[1,j]) and (alt ge dat[1,j+1]))) do j = j + 1
    tval = dat[0,j] + (dat[0,j+1] - dat[0,j]) / (dat[1,j+1] - dat[1,j]) * (dat[1,j+1] - alt)
    ttick = [!x.crange[0], tval, !x.crange[1]]
    if alt mod 50.0 then begin
      tlen = 0.01
      if alt ge 30 then begin
        tstring = ' '
      endif else begin
        tstring = strcompress(string(alt), /remove_all)
      endelse
    endif else begin
      tlen = 0.02
      tstring = strcompress(string(alt), /remove_all)
    endelse
    if (tval gt !x.crange[0]) and (tval lt !x.crange[1]) then $
      axis, xaxis = 1, $
        xtitle = 'Altitude [km]', $
        xstyle = 1, $
        xrange = !x.crange, $
        xticks = 2, $
        xtickv = ttick, $
        xticklen = tlen, $
        xtickname = [' ', tstring, ' ']
  endif
endfor

lnfl = 0
openr, 1, 'flasar.dat'
while not eof(1) do begin
  readf, 1, dummy
  lnfl = lnfl + 1
endwhile
close, 1

datfl = dblarr(2, lnfl)

openr, 1, 'flasar.dat'
readf, 1, datfl
close, 1

tval = dblarr(lnfl)

for i = 0, lnfl-1 do begin
  if (datfl[0,i] le dat[1,0]) and (datfl[0,i] ge dat[1,ln-1]) then begin
    j = long(0.)
    while (not ((datfl[0,i] le dat[1,j]) and (datfl[0,i] ge dat[1,j+1]))) do j = j + 1
    tval[i] = dat[0,j] + (dat[0,j+1] - dat[0,j]) / (dat[1,j+1] - dat[1,j]) * (dat[1,j+1] - datfl[0,i])
  endif else tval[i] = -1
endfor

datfl[1,*] = zwd * datfl[1,*] * cos(9. * !DPI / 180.) / cos(45. * !DPI / 180.)

oplot, tval, datfl[1,*], $
  linestyle = 2
oplot, [tval[0], !x.crange[1]], [0.,0.], linestyle = 2

zws = interpol(datfl[1,*], tval, dat[0,*])

for i = long(0), ln - 1 do begin
  if dat[1,i] le 0. then zws[i] = 0.
endfor

!p.region = [0., 0., 1.0, 0.5]
plot, dat[0,*], zws - dat[15,*], $
  position = [0.1, 0., 1.0, 0.5], $
  xtitle = '!18t - t!D0!N!X [min]', $
;  ytitle = 'Residuals of Zonal Wind Speed Recovery [m/s]', $
  ytitle = 'Residuals [m/s]', $
  ystyle = 1, $
;  yrange = [-3., 3.], $
  yrange = yrres, $
;  xrange = [0., 180.], $
  xstyle = 9

for alt = 160, 0, -10 do begin
  if (alt le dat[1,0]) and (alt ge dat[1,ln-1]) then begin
    j = long(0.)
    while (not ((alt le dat[1,j]) and (alt ge dat[1,j+1]))) do j = j + 1
    tval = dat[0,j] + (dat[0,j+1] - dat[0,j]) / (dat[1,j+1] - dat[1,j]) * (dat[1,j+1] - alt)
    ttick = [!x.crange[0], tval, !x.crange[1]]
    if alt mod 50.0 then begin
      tlen = 0.01
      if alt ge 30 then begin
        tstring = ' '
      endif else begin
        tstring = strcompress(string(alt), /remove_all)
      endelse
    endif else begin
      tlen = 0.02
      tstring = strcompress(string(alt), /remove_all)
    endelse
    if (tval gt !x.crange[0]) and (tval lt !x.crange[1]) then $
      axis, xaxis = 1, $
;        xtitle = 'Altitude [km]', $
        xtitle = ' ', $
        xstyle = 1, $
        xrange = !x.crange, $
        xticks = 2, $
        xtickv = ttick, $
        xticklen = tlen, $
;        xtickname = [' ', tstring, ' ']
        xtickname = [' ', ' ', ' ']
  endif
endfor

oplot, !x.crange, [0.,0.], linestyle = 2

!p.multi = 0
!p.region = 0.
;plot, dat[0,*], zws - dat[15,*], $
;  xtitle = '!18t - t!D0!N!X [min]', $
;;  ytitle = 'Residuals of Zonal Wind Speed Recovery [m/s]', $
;  ytitle = 'Residuals [m/s]', $
;  ystyle = 1, $
;  xrange = [0., 1.], $
;  xstyle = 1
;oplot, !x.crange, [0.,0.], linestyle = 2

if (medium eq 0 or medium eq 1) then begin
  device, /close  ; close 'output.ps'
  print, 'Rename "output.ps" or "output.eps" before you call this program again !!'
endif

end