;+
; NAME:
;  rt_solread
; PURPOSE: (one line)
;  read solar spectrum and cross sections
; DESCRIPTION:
;  read solar spectrum and cross sections
; CATEGORY:
;  RT
; CALLING SEQUENCE:
;  rt_solread, wave, flux, crs, nsol
; INPUTS:
;  None
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
;  cycle = 0 for solar max, 1 for solar min
;  refdir: directory of soluvmax.dat or soluvmin.dat
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  wave[nsol] - wavelength in A
;  dwave[nsol] - width of wavelength interval in A
;  flux[nsol] - flux in photon/cm^2/s/A
;  crs[nsol,4] - cross section of N2, CH4, C2H2, C2H6 (cm^2)  
;  nsol - number of solar lines
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2004 Mar 20, Leslie Young SwRI
;  based on solread, written by Roger Yelle, 1990
;  2004 Mar 21 LAY added dwave, changed flux to 
;   ph/cm^2/s/A, not ph/cm^2/s/interval
;-
pro rt_solread, wave, dwave, flux, crs, nsol, cycle = cycle, refdir = refdir


;=================
; set up filename
if not keyword_set(refdir) then begin
    if !version.os_family eq 'unix' then begin
        refdir = '~/reference/solspec/solspecUV_rvy/'
    endif else begin
        refdir = 'Laputa:Users:layoung:reference:solspec:solspecUV_rvy:'
    end
end
if not keyword_set(cycle) then cycle = 0 
if cycle eq  0 then fn = 'soluvmax.dat' else fn = 'soluvmin.dat'
fn = refdir+fn

;=================
; read it in
readcol, fn, wave, flux, crs0, crs1, crs2, crs3, skip=1, /silent
crs = [ [crs0], [crs1], [crs2], [crs3] ]
nsol = n_elements(wave)

dwave = (wave[2:nsol-1]-wave[0:nsol-3])/2
dwave = [dwave[0], dwave, dwave[nsol-3]]
flux = flux / dwave

end

pro rt_solreadTEST
  print, 'hi'
  rt_solread, wave, dwave, flux, crs, nsol, cycle=1

;  readcol,$
;    'Laputa:Users:layoung:reference:solspec:solspecUV_gladstone:solspec.dat', $
;    wave2, flux2
;  readcol,$
;    'Laputa:Users:layoung:reference:solspec:solarspectra:newsol.dat', $
;    wave3, flux3, eflux3

  !p.multi=[0,0,2]
  plot, wave, flux, tit='flux',xtit='wavelength (A)', $
    ytit='flux (phot/cm^2/s/A)', /ylog, ps=4
;  oplot, wave2, flux2, li=2
;  oplot, wave3, flux3, li=2
  plot, wave, crs[*,0], tit='',xtit='wavelength (A)', $
    ytit='cross section (cm^2)', yr=[0,max(crs)]
  for i=1,3 do  oplot, wave, crs[*,i], li=i
end