PRO sac_peak_atlas, w0, w1, wavearr, flux, intflux, FILE=file,$
                    THEKAEKARA=thekaekara,$
                    INTEGRATED = integrated, $
                    PLOTIT=plotit, EPS=eps,$
                    WM2NM=wm2nm


;Reads the ASCII file with the data from the Beckers, Bridges, and
;Gilliam 1976 Solar Flux Atlas. Printed version is AFGL-TR-76-0126(I):
;"A High Resolution Spectral Atlas of the Solar Irradiance from 380 to
;700 Nanometers", June 1976.

;ASCII file was obtained from R. Kurucz's website on
;27-Oct-2005: http://kurucz.harvard.edu/sun/sacpeak/

;INPUTS:
; W0 = starting wavelength IN NANOMETERS in format 430.5
; W1 = ending wavelength IN NANOMETERS in same format.
;
; Granularity is 0.1 nm (1 Angstrom). Within each 1A bandpass,
; there are 200 flux numbers in the atlas giving a spectral resolution of 0.01 A.

;KEYWORDS:
; THEKAEKARA: if set, returns the flux matched to Thekaekara's
; spectrum in Thekaekara, M.P., Applied Optics, 13, 518, 1974.
; INTEGRATED: returns the integrated flux over the wavelength range.
; PLOTIT: if set, plot the results.
; EPS: if set, produce an EPS file of the plot.
; WM2NM: if set, return flux in W/m^2/nm units.
; 
; Default behavior: return the flux matched to Labs & Neckel's
; spectrum in Labs, D., and Neckel, H., Solar Physics, 15, 79, 1970.

;OUTPUTS
; WAVEARR = Floating point array of wavelength values for each 0.1 nm group.
; FLUX = Floating point array of flux in Watts/m^2/micron matched to
; either Labs&Neckel or Thekaekara. If keyword WM2NM is set, units are
; Watts/m^2/nm instead.

;Data file format:
;200 integer flux values in a 16x12 array with a final row of 8 flux
;values followed by a number of the form "38002533919794". This number
;is the starting wavelength in A (3800) followed by the Thekaekara
;conversion factor (2.5339) followed by the Labs&Neckel conversion
;factor (1.9794). Ugh - who the hell came up with this lovely format????

;NOTE: although the Labs & Neckel and Thekaekara factors are supposed
;to convert the flux values to those observed outside the
;Earth's atmosphere, the spectrum still includes terrestrial
;absorption lines. Therefore use caution when calculating integrated
;flux over a wavelength interval - the terrestrial lines will decrease
;this value relative to a true space-based spectral irradiance
;measurement.

if KEYWORD_SET(file) then afile = file else $
  afile = '$SSW/hinode/sot/idl/fg/util/sacpeak.asc'

if w0 lt 380 then begin
    MESSAGE,'Minimum Sac Peak Atlas wavelength = 380 nm. Please enter w0 >= 380 nm.', /INFORM
    RETURN
end

if w1 gt 699.9 then begin
    MESSAGE,'Maximum Sac Peak Atlas wavelength = 699.9 nm. Please enter w1 <= 699.9 nm.', /INFORM
    RETURN
end


f16 = INTARR(16)
wstring = ' '

GET_LUN,lun
PRINT,'Opening ',afile
OPENR,lun,afile

flux = -1
;wavearr = -1.0
wave = 380.0
;Read the file until we find the starting wavelength
while w0 gt wave do begin
    f200 = -1
    for j=1,12 do begin
        READF,lun,f16
        f200 = [f200,f16]
    end
    READF,lun,wstring
    s = STRSPLIT(wstring,' ',/extract)
    for i=0,7 do f200 = [f200,FIX(s[i])]
    wave = FLOAT(STRMID(s[8],0,4))/10.   ;wavelength in nm
    thfac = FLOAT(STRMID(s[8],4,5))/10000.
    lnfac = FLOAT(STRMID(s[8],9,5))/10000.
    if KEYWORD_SET(thekaekara) then f200 *= thfac else f200 *= lnfac
end
;Now have 1st block of data in f200
;wavearr = [wavearr,wave]
flux = [flux,f200[1:*]]   ;get rid of leading -1 in f200

while wave lt w1 do begin
    f200 = -1
    for j=1,12 do begin
        READF,lun,f16
        f200 = [f200,f16]
    end
    READF,lun,wstring
    s = STRSPLIT(wstring,' ',/extract)
    for i=0,7 do f200 = [f200,FIX(s[i])]
    wave = FLOAT(STRMID(s[8],0,4))/10.   ;wavelength in nm
    thfac = FLOAT(STRMID(s[8],4,5))/10000.
    lnfac = FLOAT(STRMID(s[8],9,5))/10000.
    if KEYWORD_SET(thekaekara) then f200 *= thfac else f200 *= lnfac
;    wavearr = [wavearr,wave]
    flux = [flux,f200[1:*]]
end

FREE_LUN,lun

;wavearr = wavearr[1:*]
flux = flux[1:*]
wavearr = FINDGEN(N_ELEMENTS(flux))*0.0005 + w0
good = WHERE(wavearr ge w0 and wavearr le w1)
flux = flux[good]
wavearr = wavearr[good]

if KEYWORD_SET(wm2nm) then flux = flux/1000.

if KEYWORD_SET(plotit) then begin
    if KEYWORD_SET(eps) then begin
        SET_PLOT,'ps'
        ps_reset
        !P.FONT=0
        DEVICE,/encap, /preview, /portrait, /helvetica, /color,bits=8    
        DEVICE,filename='sacpeak_atlas.eps'
        DEVICE,/helvetica
    end
    if KEYWORD_SET(wm2nm) then yt = 'Flux, W/m^2/nm' else yt='Flux, W/m^2/micron'
    PLOT,wavearr,flux,psym=3,$
         xtitle='Wavelength, nm',ytitle=yt,$
         title='Solar Spectral Irradiance. NSO Sac Peak Atlas, 1976'
    OPLOT,wavearr,flux
    if KEYWORD_SET(eps)then begin
        DEVICE,/close
        SET_PLOT,'x'
    end
end

if KEYWORD_SET(integrated) then begin
;Integrated flux calculation. Flux is in units of Watts/m^2/micron
;This is the same number value as flux in erg/sec/cm^2/nm.
    if N_ELEMENTS(wavearr) gt 10 then intflux = INT_TABULATED(wavearr,flux) else $
      intflux = TOTAL(flux*0.0005)
end

RETURN
END