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='rho_n.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='rho_n.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

xyouts, 0.5, 0.5, /normal, '!5'

dummy = ' '

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

dat_iono = dblarr(3,401)
openr, 1, 'n_iono.dat'
readf, 1, dat_iono
close, 1

dat_neut = dblarr(3,401)
openr, 1, 'n_neut.dat'
readf, 1, dat_neut
close, 1

ymax = max([max(dat_iono[1,*]), max(dat_neut[1,*])])
ymin = min([min(dat_iono[1,*]), min(dat_neut[1,*])])

;plot, dat_iono[0,*], dat_iono[1,*], /ylog, $
plot, dat_iono[1,*], dat_iono[0,*], /xlog, $
  ytitle = 'Altitude [km]', $
  xtitle = 'Density [m!U-3!N]', $
;  yrange = [0., 200.], $
  ystyle = 1, $
  xrange = [ymin, ymax], $
  xstyle = 3, $
  linestyle = 1

;oplot, dat_neut[0,*], dat_neut[1,*]
oplot, dat_neut[1,*], dat_neut[0,*]

dat_iono[2,*] = -1. * dat_iono[2,*]

ymax = max([max(dat_iono[2,*]), max(dat_neut[2,*])])
ymin = min([min(dat_iono[2,*]), min(dat_neut[2,*])])

;plot, dat_iono[0,*], dat_iono[2,*], /ylog, $
plot, dat_iono[2,*], dat_iono[0,*], /xlog, $
  ytitle = 'Altitude [km]', $
  xtitle = 'Refractivity', $
;  yrange = [0., 200.], $
  ystyle = 1, $
  xrange = [ymin, ymax], $
  xstyle = 3, $
  linestyle = 1

;oplot, dat_neut[0,*], dat_neut[2,*]
oplot, dat_neut[2,*], dat_neut[0,*]

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

!x.thick = 1.
!y.thick = 1.
!p.charsize = 1.
!p.charthick = 1.

end