; Takes Bus-DeMeo taxonomy individual spectra and converts them to VISTA colors
; Returns flux and errors


function spec2vista, specFile, vistaDIR, plotTAXO=plotTAXO, dirFig=dirFig



  xRang = [0.3, 2.5]
  yRang = [0.5, 2.1]


;--------------------------------------------------------------------;
;--------------------------------------------------------------------;
;-- TAG --  1 --   Initialization, Constants...                  ----;
;--------------------------------------------------------------------;
;--------------------------------------------------------------------;  



  readcol, specFile, ws, fs, es, format='(F,F,F)',/silent
    g = where( fs ne -1 )
  ws=ws(g)
  fs=fs(g)
  es=es(g)

  ord=sort(ws)
  ws=ws(ord)
  fs=fs(ord)
  es=es(ord)

 ;-                                               Z      Y   J       H       K
  if not keyword_set(filtwave) then filtwave = [0.877,1.020,1.252,1.645,2.147]


  flux   = fltarr(5)
  fluxerr= fltarr(5)

  vistaTAXO = {Z: 0., $
               Y: 0., $
               J: 0., $
               H: 0., $
               K: 0. }


 ;-Read VISTA QE
  readcol, vistaDIR+'qe.tab', wQE, tQE, /silent
  readcol, vistaDIR+'VISTA_M1_Reflectivity_forETC_2009-Sep12.txt', wM1, tM1, /silent
  readcol, vistaDIR+'VISTA_M2_Reflectivity_forETC_2007-Jun19.txt', wM2, tM2, /silent


 ;-Read VISTA Filter Response
  readcol, vistaDIR+'VISTA_Filters_at80K_forETC_Z.dat',  wZ, tZ, /SILENT
  readcol, vistaDIR+'VISTA_Filters_at80K_forETC_Y.dat',  wY, tY, /SILENT
  readcol, vistaDIR+'VISTA_Filters_at80K_forETC_J.dat',  wJ, tJ, /SILENT
  readcol, vistaDIR+'VISTA_Filters_at80K_forETC_H.dat',  wH, tH, /SILENT
  readcol, vistaDIR+'VISTA_Filters_at80K_forETC_Ks.dat', wK, tK, /SILENT



 ;-WAvelength convertion to microns
  wZ /= 1.e3 &  wY /= 1.e3 &  wJ /= 1.e3 &  wH /= 1.e3 &  wK /= 1.e3
  wQE /= 1.e3  & wM1 /= 1.e3 &  wM2 /= 1.e3

  tZ( where( wZ le 0.78 or wZ ge 1.00) ) = 0
  tY( where( wY le 0.80 or wY ge 1.20) ) = 0
  tJ( where( wJ le 1.05 or wJ ge 1.50) ) = 0
  tH( where( wH le 1.48 or wH ge 2.00) ) = 0
  tK( where( wK le 1.88 or wK ge 2.45) ) = 0

  validZ = where( tZ ne 0 )
  validY = where( tY ne 0 )
  validJ = where( tJ ne 0 )
  validH = where( tH ne 0 )
  validK = where( tK ne 0 )




;plot, ws, fs, yr=[0,2]
;oplot, wZ, tZ/100.
;oplot, wY, tY/100.
;oplot, wJ, tJ/100.
;oplot, wH, tH/100.
;oplot, wK, tK/100.
;stop


;--------------------------------------------------------------------;
;--------------------------------------------------------------------;
;-- TAG --  2 --   Compute Fluxes, Colors...                     ----;
;--------------------------------------------------------------------;
;--------------------------------------------------------------------;  

  normaI = 2
  cutTrans = 25.

  specZ = interpol( fs, ws, wZ )
  specY = interpol( fs, ws, wY )
  specJ = interpol( fs, ws, wJ )
  specH = interpol( fs, ws, wH )
  specK = interpol( fs, ws, wK )

  errZ = interpol( es, ws, wZ )
  errY = interpol( es, ws, wY )
  errJ = interpol( es, ws, wJ )
  errH = interpol( es, ws, wH )
  errK = interpol( es, ws, wK )

  CASE normaI OF
     0: BEGIN
        normaR = where( tZ ge cutTrans )
        normaV = mean(specZ( normaR ))
     END
     1: BEGIN
        normaR = where( tY ge cutTrans )
        normaV = mean(specY( normaR ))
    END
     2: BEGIN
        normaR = where( tJ ge cutTrans )
        normaV = mean(specJ( normaR ))
     END
     3: BEGIN
        normaR = where( tH ge cutTrans )
        normaV = mean(specH( normaR ))
     END
     4: BEGIN
        normaR = where( tK ge cutTrans )
        normaV = mean(specK( normaR ))
     END


  ENDCASE

  fs /= normaV
  es /= normaV

  specZ /= normaV
  specY /= normaV
  specJ /= normaV
  specH /= normaV
  specK /= normaV


 ;--2.1-- Convert DeMeo Spectra into VISTA Colors (Normalized to R)
  fluxes=[ total( specZ(validZ)*tZ(validZ) ) / total(tZ(validZ)), $
           total( specY(validY)*tY(validY) ) / total(tY(validY)), $
           total( specJ(validJ)*tJ(validJ) ) / total(tJ(validJ)), $
           total( specH(validH)*tH(validH) ) / total(tH(validH)), $
           total( specK(validK)*tK(validK) ) / total(tK(validK))  ]
  colors = -2.5*alog10( fluxes/fluxes(normaI) )
    vistaTAXO.Z = colors(0)
    vistaTAXO.Y = colors(1)
    vistaTAXO.J = colors(2)
    vistaTAXO.H = colors(3)
    vistaTAXO.K = colors(4)

  flux(*) = [ 10.^( -0.4* (vistaTAXO.Z) ), $
              10.^( -0.4* (vistaTAXO.Y) ), $
              10.^( -0.4* (vistaTAXO.J) ), $
              10.^( -0.4* (vistaTAXO.H) ), $
              10.^( -0.4* (vistaTAXO.K) ) ]
    

 ;--2.2-- Associate DeMeo Deviation into VISTA Colors
  fluxerr(*)=[ mean(errZ(validZ)), $
               mean(errY(validY)), $
               mean(errJ(validJ)), $
               mean(errH(validH)), $
               mean(errK(validK))  ]




  if keyword_set(plotTAXO) eq 1 then begin 

   ;---------------------------------------------------------------------------
   ;--2.3-- TaxoPlot
   ;--2.3.1-- Preparation
    set_plot,'ps'
    device, xs=12, ys=8

    split = strsplit( specFile,'/',/extract, count=NN)
    nameFig = strmid( split(NN-1),0,strlen(split(NN-1))-3)+'eps'

    device, filename=dirFig+nameFig, /color, encapsulated=0
    !P.CHARSIZE = 0.8
    !X.STYLE=1 & !Y.STYLE=1
    !Y.TICKINTERVAL=.5 &   !X.TICKINTERVAL=.2
    !Y.TICKLEN=.03
    !X.TICKLEN=.05


    plot, ws, fs, xr=xRang, yr=yRang, xst=1,yst=1,/NODATA, $
          xtickint=0.5, xminor=5


   ;--2.3.2-- Taxonomic class range
    range=where( ws ge 0.35 and ws le 2.4 )
    xClass = [ reform(ws(range)), $
               reform(reverse(ws(range)))]
    yClass = [ reform(fs(range)-es(range)), $
               reform(fs(reverse(range))+es(reverse(range))) ]
  ;  polyfill, xClass, yClass, /DATA, color=200

    oplot, ws, fs
    oplot, xClass, yClass, color=130

   ;--2.3.3-- VISTA Color for Taxonomy Classes
;    define_psym, 1, symbol=11
;    oploterror, filtwave, flux(*), fluxerr(*), psym=3
    cgPlot, /over, filtwave, flux(*), psym='filled circle'




   ;--2.3.4-- Class Name
;    xyouts,0.5,1.5, demeoTAXO.class(kclass), charsize=1.5, charthick=2

    device, /close
    set_plot, 'X'
  endif


  
  return,  {ref: flux, $
            unc: fluxerr, $
            col: colors}
end