;+
; NAME:
;   desi_quicklya
;
; PURPOSE:
;   Compute the Lyman-alpha S/N for the mean QSO at different g-mag, z
;
; CALLING SEQUENCE:
;   desi_quicklya
;
; INPUTS:
;
; OPTIONAL INPUTS:
;   exptime    - Exposure time; default to 4000 sec
;
; OUTPUTS:
;
; OPTIONAL OUTPUTS:
;
; COMMENTS:
;   Input file has wavelength and mean QSO spectra at z=2.0-4.75 spaced
;   every 0.25 in redshift:
;     $DESIMODEL/data/spectra/spec-lya.dat
;   Write one output files per magnitude slice, each containing the
;   S/N per Ang for all redshifts:
;     $DESIMODEL/data/spectra/sn-spec-lya-XX.XX.dat
;
; EXAMPLES:
;
; BUGS:
;
; DATA FILES:
;
; REVISION HISTORY:
;   01-Jul-2014  Written by D. Schlegel, LBL
;-
;------------------------------------------------------------------------------
pro desi_quicklya, exptime=exptime1

   if (keyword_set(exptime1)) then exptime = exptime1 $
    else exptime = 4000.

   topdir = getenv('DESIMODEL')
   if (NOT keyword_set(topdir)) then $
    message, 'Enviroment DESIMODEL must be set'

   ; Read in mean Ly-a spectra at z=2.0 - 4.75 in steps of 0.25
   infile = topdir+'/data/spectra/spec-lya.dat'
   readcol, infile,wave,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12
   flux = [[f1],[f2],[f3],[f4],[f5],[f6],[f7],[f8],[f9],[f10],[f11],[f12]]
   nobj = 12

   ; Step through different normalizations of the magnitude
   for mag=19.25, 22.75, 0.50 do begin
      for iobj=0, nobj-1 do begin

         print, 'Working on MAG=', mag, ' spectrum# ', iobj

         wavevec=10d0^(3.5d0+dindgen(6000)*1e-4)
         tmpflux=interpol(flux[*,iobj],wave,wavevec)
         flambda2fnu = wavevec^2 / 2.99792e18
         fthru = filter_thru(tmpflux * flambda2fnu, waveimg=wavevec, /toair)
         spectroflux = transpose(fthru) * 10^((22.5 + 48.6 - 2.5*17.)/2.5)
         thismag = 22.5-2.5*alog10(spectroflux)

         ; Normalize to g-band magnitude MAG
         newflux = flux[*,iobj] * 10^((thismag[1] - mag)/2.5)

         desi_quicksim, model='star', wave, newflux, exptime=exptime, $
          simdat=simdat

         ; Bin the wavelengths to 1 Ang and compute the S/N per 1 Ang
         npix = n_elements(simdat.wavelength)
         wave1 = rebin(simdat.wavelength, npix/2)
         snvec1 = rebin(simdat.snvec, npix/2) * sqrt(2)

         if (iobj EQ 0) then snvec = fltarr(npix/2, nobj)
         snvec[*,iobj] = snvec1

      endfor
      magstr = string(mag,format='(f5.2)')
      outfile = topdir+'/data/spectra/sn-spec-lya-'+magstr+'.dat'
      openw, olun, outfile, /get_lun
      printf, olun, '# Lyman-alpha forest S/N per Ang for mean quasar with mean forest'
      printf, olun, '# INFILE= '+infile
      printf, olun, '# GMAG= '+strtrim(mag,2)
      printf, olun, '# EXPTIME= '+strtrim(exptime,2)
      printf, olun, '#'
      printf, olun, '# Wave SN(z=2.0) SN(z=2.25) SN(z=2.50) SN(z=2.75) SN(z=3.00) SN(z=3.25) SN(z=3.50) SN(z=3.75) SN(z=4.00) SN(z=4.25) SN(z=4.50) SN(z=4.75)'
      for i=0L, n_elements(wave1)-1L do $
       printf, olun, wave1[i], snvec[i,*], format='(f8.2,12f6.2)'
   endfor

   return
end
;------------------------------------------------------------------------------