;
function getbins,file
;
; Read given file to get Guy Brassier wavelength bins
;
openr,lu,file,/get_lun
hdr1 = ''
readf,lu,hdr1
;
; Get number of bins:
;
nbins = 0
while (not EOF(lu)) do begin
  readf,lu,nbins
endwhile
;
; Bin structure to be returned:
;
getbins = {binstruct,         $
  ibin     :intarr(nbins),    $
  wv_low   :fltarr(nbins),    $
  wv_high  :fltarr(nbins),    $
  wv_del   :fltarr(nbins)     }
;
point_lun,lu,0
readf,lu,hdr1
ibin = 0
print,'Getbins: reading wavelength bins from file ',file
for i = 0,nbins-1 do begin
  readf,lu,ibin,wv_low,wv_high,wv_del
  getbins.ibin(i) = ibin
  getbins.wv_low(i) = wv_low
  getbins.wv_high(i) = wv_high
  getbins.wv_del(i) = wv_del
; print,format="(i5,2f12.2,f9.2)",getbins.ibin(i),getbins.wv_low(i),$
;   getbins.wv_high(i),getbins.wv_del(i)
endfor
free_lun,lu
return,getbins
end
;-----------------------------------------------------------------------
function getfluxday, file,chsolar,datfile=datfile,ps=psfile
;
; Read a Solstice ascii file containing high resolution
; spectra (.05 nm "hires" files) for a single uars day, 
; and bin into Buy Brassseur's 169 wavelength bins.
; (Solstice file obtained from anony ftp rescha.colorado.edu,
;  /pub/solstice)
;
; E.g.: getflux,"sol_hires_0582.dat",582
;
print,' '
print,'Getflux: data file=',file
if keyword_set(datfile) gt 0 then print,'datfile=',datfile
if keyword_set(ps) then print,'Will make ps file'
;
; Read_dat.pro is from the same ftp site as the data, see above.
;
uars_day = -1
data = read_dat(file,day=uars_day)
date = uday_to_date(uars_day)
print,'UARS Day = ',uars_day,' Date = ',date
;
wavestep = 0.05
factor = 503.4020               ; = 1.E-13 / h / c
photon_flux = data(1,*) * factor * data(0,*)
;
; Get Brasseur wavelength bins (see getbins function above)
;
bins = getbins("brasseur_bins.dat")
nbins = n_elements(bins.ibin)
flux = fltarr(nbins)
count = 0
;
;plot,bins.wv_low
;oplot,bins.wv_high
;if (n_elements(psfile) eq 0) then cursor,x,y,/wait
;
; Loop through Brasseur bins, summming fluxes in each
; (convert flux to photons/cm2/s)
;
for ibin = 0,nbins-1 do begin
  ww = where((data(0,*) ge bins.wv_low(ibin)) and $
             (data(0,*) lt bins.wv_high(ibin)),count)
  if (count gt 0) then begin
    flux(ibin) = total(photon_flux(ww)) * wavestep
  endif else begin
    flux(ibin) = 0.
  endelse
; print,format="('ibin=',i4,'  low,high=',2f8.2,'  flux=',e12.4)",$
;   ibin,bins.wv_low(ibin),bins.wv_high(ibin),flux(ibin)
endfor
iflux = where(flux ne 0.)
pltflux = flux(iflux)
;
; Calculate min,max,mean: 
;
;fluxmin = min(flux(where(flux gt 0.)),max=fluxmax)
fluxmin = min(pltflux,max=fluxmax)
fluxave = total(pltflux)/float(n_elements(iflux))
print,'Flux min,max,mean = ',fluxmin,fluxmax,fluxave
title1 = "UARS DAY " + string(format="(i5)",uars_day) + "  " + date + $
         " " + strtrim(chsolar,2)
title1 = strtrim(title1,2)
title2 = " Flux min,max,mean = " + string(format="(3e9.2)",$
         fluxmin,fluxmax,fluxave)
title2 = strtrim(title2,2)
;
; Write to ascii file if requested:
;
if keyword_set(datfile) gt 0 then begin
  openw,luwr,datfile,/get_lun
  printf,luwr,title1
  printf,luwr,format="('  bin      wv_low      wv_high    wv_del   flux')"
  format="(i5,2f12.2,f9.2,e15.4)"
  for i=0,nbins-1 do begin
    printf,luwr,i+1,bins.wv_low(i),bins.wv_high(i),$
      bins.wv_del(i),flux(i),format=format
  endfor
  close,luwr
  print,'Wrote file ',datfile
  free_lun,luwr
endif
;
; Open ps file if requested:
;
if (n_elements(psfile) gt 0) then begin
  set_plot,'PS'
  device,filename=psfile
endif else begin
  set_plot,'X'
endelse
;
; Plot:
;
ibins = indgen(nbins)
ibins = ibins(iflux)+1
log = 1 ; plot log10 of flux
binmin = min(bins.wv_low(iflux),max=binmax)
xtitle = 'BRASSEUR BINS (min,max wavelengths = ' + $
  strtrim(string(binmin,format="(f7.2)"),2) + ', ' + $
  strtrim(string(binmax,format="(f7.2)"),2) + ' nm)'
position = [.15,.25,.88,.88]
if (log gt 0) then begin
  plot_io,ibins,pltflux,xtitle=xtitle,position=position,/norm,charsize=1.3,$
    ytitle='SOLSTICE FLUX (photons/s/cm^2)',xmargin=[13,5]
endif else begin
  plot,ibins,pltflux,xtitle=xtitle,position=position,/norm,charsize=1.3,$
    ytitle='SOLSTICE FLUX (photons/s/cm^2)',xmargin=[13,5]
endelse
;
; Crude extra x-axis showing selected wavelenths (nm):
;
psfac = 0.
if (n_elements(psfile) gt 0) then psfac = .04
ypos = !y.window(0)-0.11-psfac
plots,!x.window(0),ypos,/norm
plots,!x.window(1),ypos,/continue,/norm
lblbins = [1,10,20,30,40,50,60,70,80,90,100,110,120]
for i=0,n_elements(lblbins)-1 do begin
  lamda = 0.5*(bins.wv_low(lblbins(i)-1)+bins.wv_high(lblbins(i)-1))
  axpos = convert_coord(lblbins(i),0.,/data,/to_norm)
  xpos = axpos(0)
;
; Tic mark:
  plots,xpos,ypos,/norm
  plots,xpos,ypos-.01,/norm,/continue
;
; Tic Label:
  label = strtrim(string(lamda,format="(f7.2)"),2)
  xyouts,xpos,ypos-.02,label,/norm,orientation=90,alignment=1.0
endfor
;
; Axis label:
xyouts,0.5*(!x.window(0)+!x.window(1)),ypos-.12-psfac,"Wavelength (nm)",$
  /norm,charsize=1.3,alignment=0.5
;
; Label selected bins:
; (e.g., Lyman-alpha (bin #8, 121.2-122.0))
;
;lblbins = [8,19,23,29,36,42]
lblbins = [8,19,23]
for i=0,n_elements(lblbins)-1 do begin
  plots,float(lblbins(i)),flux(lblbins(i)-1),psym=4
  locnorm = convert_coord(float(lblbins(i)),flux(lblbins(i)-1),$
    /data,/to_norm)
  lamda = 0.5*(bins.wv_low(lblbins(i)-1)+bins.wv_high(lblbins(i)-1))
  label = strtrim(string(flux(lblbins(i)-1),format="(e10.2)"),2) + " (" + $
          strtrim(string(lamda,format="(f7.2)"),2) + " nm)"
  xyouts,locnorm(0),locnorm(1)+.015,label,/norm,orientation=90
endfor
;
; Two-part title at top:
;
xpos = 0.5*(!x.window(0)+!x.window(1))
ypos = !y.window(1)+.07
xyouts,xpos,ypos,title1,alignment=0.5,/norm,charsize=1.3
ypos = !y.window(1)+.03
xyouts,xpos,ypos,title2,alignment=0.5,/norm,charsize=1.3
if (n_elements(psfile) eq 0) then cursor,x,y,/wait
;
; Close ps file:
;
if (n_elements(psfile) gt 0) then begin
  device,/close_file
  print,'Wrote plot to file ',psfile
endif
;
; Return flux:
;
return,flux
end
;-----------------------------------------------------------------------
pro getflux
;
; Call getfluxday for solar minimum, medium, and maximum uars days:
;
; At the ftp site ftp rescha.colorado.edu (/pub/solstice), there are
;   ascii "hires" files for solar cycle 22:
;
; sol_hires_0144.dat: 
;   Reference solar UV irradiance spectrum for solar cycle 22 
;   maximum value.
; sol_hires_0200.dat:
;   Reference solar UV irradiance spectrum for solar cycle 22
;   average maximum condition.
; sol_hires_0582.dat:
;   Reference solar UV irradiance spectrum for solar cycle 22
;   moderate condition.
; sol_hires_1209.dat:
;   Reference solar UV irradiance spectrum for solar cycle 22
;   average minimum condition. 
;
; Calls for X-display with no output data files:
;
;flux = getfluxday("sol_hires_1209.dat","(Solar Minimum)")
;flux = getfluxday("sol_hires_0582.dat","(Solar Medium)")
;flux = getfluxday("sol_hires_0144.dat","(Solar Maximum)")
;
; File for fortran data statements:
;
datafile = "smin_med_max.data"
openw,luout,datafile,/get_lun
;
; Solar minimum:
;
flux = getfluxday("sol_hires_1209.dat","(Solar Minimum)", $
  datfile="smin_day1209.dat",ps="smin_day1209.ps")
wrdata,flux,"smin",luout=luout
;
; Solar medium:
;
flux = getfluxday("sol_hires_0582.dat","(Solar Medium)",  $
  datfile="smed_day0582.dat",ps="smed_day0582.ps")
wrdata,flux,"smed",luout=luout
;
; Solar maximum:
;
flux = getfluxday("sol_hires_0144.dat","(Solar Maximum)", $
  datfile="smax_day0144.dat",ps="smax_day0144.ps")
wrdata,flux,"smax",luout=luout
;
; Make data statements for time-gcm:
;
print,'wrdata wrote data statements file ',datafile
free_lun,luout
;
end