;-- same soft as euclid-demeo, but limited to 11 classes
; O R are dropped
; complexes are merged together



;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--I.1-- Directories
  spawn, 'hostname', host
  case strTrim(host) of
     'hyperion': dirEuclid = '/observ/euclid/'
     'endymion': dirEuclid = '/home/bcarry/work/data/euclid/'
     else: stop
  endcase
  filterDir = dirEuclid+'filters/'


  kboDIR = '/astrodata/Catalog/Taxonomy/2016-Merlin/'

  ecCol  = ['Black', 'Orange', 'Cornflower blue', 'Firebrick']

  doTABLE=1


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Filter/Spectra to Colors                           -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--II.1-- Euclid filters
  filtName=['VIS', 'Y', 'J', 'H']
  filtWave=[0.725,1.020,1.252,1.645]
  nbFilt=n_elements(filtName)
  
  centWave = fltarr(4)
  waveToMicron=1e-4


 ;--II.2-- Typical Filter
  readcol, filterDIR+'EuclidVisFilterMin.txt', w, t, /SILENT
  ref = {wave:w*waveToMicron, trans:t}

 ;--II.3-- Read Bus-DeMeo average spectra
  demeo  = readBusDeMeo()
  sel = [0,1,2,7,8,9,11,13,18,19,20]
  nbSel = n_elements(sel)

  newDM = { N:nbSel, wave: demeo.wave, $
           class: replicate( {name:'', spec:demeo.class[0].spec, dev:demeo.class[0].spec}, nbSel ) }
  for k=0, nbSel-1 do begin
    newDM.class[k].name = demeo.class[sel[k]].name
    newDM.class[k].spec = demeo.class[sel[k]].spec
    newDM.class[k].dev  = demeo.class[sel[k]].dev 
  endfor
  demeo=newDM





 ;--II.4-- Read Merlin KBO average spectra
  kboName = ['BB','BR','IR','RR']
  nbKBO = 4

  readcol, kboDir+kboName[0]+'.csv', kw, kr, kd, delimiter=',', format='(F,F,F)',/Silent
  kbo = { N:nbKBO, wave: kw, $
          class: replicate( {name:'', spec:kr, dev:kd}, nbKBO ) }

  for k=0, nbKBO-1 do begin
     readcol, kboDir+kboName[k]+'.csv', kw, kr, kd, delimiter=',', format='(F,F,F)',/Silent
     kbo.class[k].name = kboName[k]
     kbo.class[k].spec = kr
     kbo.class[k].dev  = kd/2.
  endfor


  nbTot = demeo.N+kbo.N



  taxoEuclid = replicate({class:'', color:replicate({ref:0., unc:0.},nbFilt)},nbTot)
  taxoEuclid[0:demeo.N-1].class     = demeo.class.name
  taxoEuclid[demeo.N:nbTot-1].class = kbo.class.name


 ;--II.4-- Convert Spectra into Colors
  print, 'Convertion to colors'
  filt=replicate(ref,nbFilt)
  for kFilt=0, nbFilt-1 do begin

   ;--II.4.1-- Read current filter
    print, filtName[kFilt]
    readcol, filterDIR+'Euclid'+filtName[kFilt]+'FilterMin.txt', w, t, /SILENT

   ;--II.4.2-- Correct bad J Filter transmission
    if strcmp(filtName[kFilt],'J') then begin
      band = where( w ge 1.2/waveToMicron and w le 1.55/waveToMicron )
      t*=0
      t[band]=0.4
    endif

   ;--II.4.3-- Convert to micron
    filt[kFilt].wave=w*waveToMicron
    filt[kFilt].trans= t
    cur = {wave:w*waveToMicron, trans:t}
    
;    forprint, w, t, format='(F,",",F)', /Silent, comment='Wavelength,Transmission',$
;              textout='/home/bcarry/data/mining/filters/euclid/euclid/euclid.'+filtName[kFilt]+'.csv'


   ;--II.4.4-- Store mid-filter wavelength
    noZero = where( cur.trans ne 0 )
    centWave[kFilt] = mean( cur.wave[noZero] ) 

   ;--II.4.5-- Convolve spectra with filter transmission curve
;    conv = convertTaxoToColor( demeo, cur, ref )
;    taxoEuclid[*].color[kFilt].ref = conv.ref
;    taxoEuclid[*].color[kFilt].unc = conv.unc


    conv = convertTaxoToColor( demeo, cur, ref )
    taxoEuclid[0:demeo.N-1].color[kFilt].ref = conv.ref
    taxoEuclid[0:demeo.N-1].color[kFilt].unc = conv.unc

    conv = convertTaxoToColor( kbo, cur, ref )
    taxoEuclid[demeo.N:nbTot-1].color[kFilt].ref = conv.ref
    taxoEuclid[demeo.N:nbTot-1].color[kFilt].unc = conv.unc



  endfor
  filtWave = centWave






;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---  III -- Multi-panel Figure                                 -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  if doTABLE eq 1 then begin

    print, 'Graphic'

   ;--III.1-- EPS/PNG output
    cgPS_open, xSize=29.7, ySize=29.7, Filename=filterDir+'euclid-taxonomy.eps', $
               /metric, /decomposed, /landscape, /encapsulated, language_level=2, /quiet

  
    xBot=0.055
    yBot=0.06
    xLen=0.225
    yLen=0.225;30
    xShi=0.00
    yShi=0.00

    nbX = 4
    nbY = 4

    cgText, /normal, 0.5,.0055, 'Wavelength ('+cgGreek('mu')+'m)', align=0.5, charsize=1.5, charthick=2
    cgText, /normal, 0.02, 0.5, 'Reflectance (normalized to VIS)',orientation=90, align=0.5, charsize=1.5, charthick=2

    xTicks=replicate(' ',10)

    xRang=[0.3,2.4]
    yRang=[0.3,2.5]



   ;--III.2-- Sub-panel: one for each class
    for kC=0, nbTot-1 do begin

;      if kC eq demeo.N then goto, j2loop


     ;--III.2.1-- Plot position
      if kC lt demeo.N then begin
        xPos = kC mod nbX
        yPos = (kC / nbX)
      endif else begin
        xPos = (kC+1) mod nbX
        yPos = ((kC+1) / nbX)
      endelse

      if xPos eq 0 then begin
        yTicks = string([0.5,1,1.5,2.],format='(F3.1)')
      endif else yTicks=replicate(' ',10)

      if yPos eq (nbY-1) then begin
        xTicks = string([0.5,1,1.5,2.,2.5],format='(F3.1)')
      endif else xTicks=replicate(' ',10)


      pPos = [xBot + xPos*(xLen+xShi), $
              yBot + (nbY-1)*(yLen+yShi)- yPos*(yLen+yShi), $
              xBot + xPos*(xLen+xShi) + xLen, $
              yBot + (nbY-1)*(yLen+yShi)- yPos*(yLen+yShi) + yLen]


      if kC lt demeo.N then begin
         wave = demeo.wave
         spec = demeo.class[kC].spec
         dev  = demeo.class[kC].dev
         name = demeo.class[kC].name+'-type'
      endif else begin
;         yRang=[0.3,3.7]
         wave = kbo.wave
         spec = kbo.class[kC-demeo.N].spec
         dev  = kbo.class[kC-demeo.N].dev
         name = kbo.class[kC-demeo.N].name+'-type'
      endelse

      norma = where( wave ge 0.55 and wave le 0.85 )
      range=where( wave ge 0.35 and wave le 2.4 )
      spec /= mean(spec[norma])

     ;--III.2.2-- Create Plot
      cgPlot, wave, spec, /NoData, /NoErase, $
              xStyle=1, xRange=xRang, xTickInt=0.5, xTickName=xTicks, xMinor=5, $
              yStyle=1, yRange=yRang, yTickInt=0.5, yTickName=yTicks, yMinor=5, $
              xTickLen=0.04, $
              yTickLen=0.03, $
              position=pPos, charSize=1.


     ;--III.2.3-- Plot DeMeo Spectra
      xClass = [ reform(wave[range]), $
                 reform(reverse(wave[range]))]
      yClass = [ reform(spec[range]-dev[range]), $
                 reform(spec[reverse(range)]+dev[reverse(range)]) ]

      cuted = where(yClass ge xRang[1])
      if cuted[0] ne -1 then yClass[cuted]=yRang[1]
      cuted = where(yClass le xRang[0])
      if cuted[0] ne -1 then yClass[cuted]=yRang[0]

      cgColorFill, xClass, yClass, color='grey'

      cgPlot, /over, wave, spec
      cgPlot, /over, xClass, yClass, color='dark grey'




     ;--III.2.4-- Overplot Photometry
;      cgErrPlot, filtWave, taxoEuclid[kC].color.ref-taxoEuclid[kC].color.unc, $
;                 taxoEuclid[kC].color.ref+taxoEuclid[kC].color.unc
;      cgPlot, /over, filtWave, taxoEuclid[kC].color.ref, psym='filled circle'
      for kFilt=0, nbFilt-1 do begin
        cgErrPlot, filtWave[kfilt], taxoEuclid[kC].color[kFilt].ref-taxoEuclid[kC].color[kFilt].unc, $
                   taxoEuclid[kC].color[kFilt].ref+taxoEuclid[kC].color[kFilt].unc, color=ecCol[kFilt], $
                   thick=3, width=.05
        cgPlot, /over, filtWave[kFilt], taxoEuclid[kC].color[kFilt].ref, psym='filled circle', color=ecCol[kFilt]
      endfor

      cgText, 0.4,2.10, name, charsize=1.5, charthick=2



j2loop:
    endfor


   ;--III.3-- Filter Transmission
    kC=demeo.N

   ;--III.3.1-- Plot position 
    xPos = kC mod nbX
    yPos = (kC / nbX)

    pPos = [xBot + xPos*(xLen+xShi), $
            yBot + (nbY-1)*(yLen+yShi)- yPos*(yLen+yShi), $
            xBot + xPos*(xLen+xShi) + xLen, $
            yBot + (nbY-1)*(yLen+yShi)- yPos*(yLen+yShi) + yLen]

   ;--III.3.2-- Create Plot
    cgPlot, wave, demeo.class[0].spec, /NoData, /NoErase, $
            xStyle=1, xRange=xRang, xTickInt=0.5, xTickName=replicate(' ',10), xMinor=5, $
            yStyle=1, yRange=[0,1], yTickInt=0.2, yTickName=yTicks, yMinor=4, $
            xTickLen=0.04, $
            yTickLen=0.03, $
            position=pPos, charSize=1.
    cgAxis, xAxis=0, yaxis=1, charSize=1.

   ;--III.3.3-- Plot filters
    for kFilt=0, nbFilt-1 do begin
      cgPlot, /over, filt[kFilt].wave, filt[kFilt].trans, color=ecCol[kFilt];'grey'
      cgText, /data, filtWave[kFilt], max(filt[kFilt].trans)+0.05, filtName[kFilt], align=0.5, color=ecCol[kFilt];'grey'
    endfor

    cgText, /data, mean(xRang),.85, 'Filter transmission', charsize=1.5, charthick=2, align=0.5





    cgPS_close, /png, Delete_PS=0

  endif


end