;--- svo service
;  wget "http://svo2.cab.inta-csic.es/theory/fps3/fps.php?Facility=HST" -O hst.xml
;  stilts tcopy in=hst.xml ifmt=votable out=hst.csv ofmt=csv
;  sed -i 's/,,/, ,/g' hst.csv

; ref: http://svo2.cab.inta-csic.es/theory/fps3/index.php?id=Paranal/FORS1.ESO1035&&mode=browse&gname=Paranal&gname2=FORS1#filter
; wget "http://svo2.cab.inta-csic.es/svo/theory/fps3/fps.php?ID=Paranal/FORS1.ESO1035" -O Paranal.FORS1.ESO1035.xml
; stilts tcopy ifmt=votable ofmt=csv in=Paranal.FORS1.ESO1035.xml out=Paranal.FORS1.ESO1035.csv


;  doSSO = 0

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

 ;--I.1-- Computer & Directories  --------------------------------------------------------------;
  spawn, 'hostname', host
  case host of
     'hyperion': root='/data/colors/WP3-Classification/'
;     'endymion': root='/home/bcarry/work/data/filters/'
;     'raul':     root='/home/bcarry/data/filters/'
     else: begin
       print, 'Where am I?'
       stop
     end
  endcase
  ssoRoot = 'sso_colors'

  ssoOut = 'bus-demeo-colors.csv'

  dirSSO = root+'SSO_spectra/'
  file_mkdir, dirSSO

;  aim='imcce'  ;- Bessell, Cousins, Johnson, SDSS, 2MASS, Gaia, Kepler
;  aim='perso'  ;- Johnson, SDSS, VISTA, Gaia, Kepler
;  aim='HST'    ;- All HST
;  aim='once'   ;- Defined below

  aim='wp3'     ;- SDSS, VISTA


 ;--I.2-- Magnitude System  --------------------------------------------------------------------;
  AB   = 0
  ST   = 0
  Vega = 01

  photSys=[keyword_set(AB), keyword_set(ST), keyword_set(Vega)]
  if total(photSys) ne 1 then begin
    print, 'Select only one Photometric System'
    stop
  endif
  if AB   eq 1 then sufSys = '-AB.csv'
  if ST   eq 1 then sufSys = '-ST.csv'
  if Vega eq 1 then sufSys = '-Vega.csv'
  ssoColor = ssoRoot+sufSys

 ;--I.3-- Definition of Instruments to Process  ------------------------------------------------;
  eSet = {facility:'', instrument:'', file:''}
  case aim of 
   ;--I.3.1-- Set for Generic Filters for VO.imcce.fr
    'imcce': begin
      nameSet='imcce'
      set=replicate(eSet,8)
      set[0].facility='Generic'
      set[0].instrument='Bessell'
      set[1].facility='Generic'
      set[1].instrument='Cousins'
      set[2].facility='Generic'
      set[2].instrument='Johnson'
      set[3].facility='SLOAN'
      set[3].instrument=' '
      set[4].facility='2MASS'
      set[4].instrument=' '
      set[5].facility='Paranal'
      set[5].instrument='VIRCAM'
      set[6].facility='GAIA'
      set[6].instrument='GAIA'
      set[7].facility='Kepler'
      set[7].instrument='Kepler'
    end
  
   ;--I.3.2-- Set for my typical filters
    'perso': begin
      nameSet='perso'
      set=replicate(eSet,5)
      set[0].facility='Generic'
      set[0].instrument='Johnson_UBVRIJHKL'
      set[1].facility='SLOAN'
      set[1].instrument=' '
      set[2].facility='Paranal'
      set[2].instrument='VIRCAM'
      set[3].facility='GAIA'
      set[3].instrument='GAIA'
      set[4].facility='Kepler'
      set[4].instrument='Kepler'
    end

   ;--I.3.3-- HST filter set
    'HST': begin
      nameSet='HST'
      set=replicate(eSet,4)
      set[0].facility='HST'
      set[0].instrument='ACS'
      set[1].facility='HST'
      set[1].instrument='NICMOS'
      set[2].facility='HST'
      set[2].instrument='WFPC2'
      set[3].facility='HST'
      set[3].instrument='WFC3'
    end

   ;--I.3.4-- SolidRock WP3 filter set
    'wp3': begin
      nameSet='wp3'
      set=replicate(eSet,2)
      set[0].facility='SLOAN'
      set[0].instrument=' '
      set[1].facility='Paranal'
      set[1].instrument='VIRCAM'
    end

   ;--I.3.4-- Local: underwork
    else: begin
      nameSet='local'
      set=replicate(eSet,4)
      set[0].facility='HST'
      set[0].instrument='ACS'
      set[1].facility='HST'
      set[1].instrument='NICMOS'
      set[2].facility='HST'
      set[2].instrument='WFPC2'
      set[3].facility='HST'
      set[3].instrument='WFC3'
    end
  endcase
  nbSet=n_elements(set)
  print, 'Number of sets: '+strTrim(string(nbSet),2)

 ;--I.X-- File Extension  ----------------------------------------------------------------------;
  sufCSV='.csv'
  sufXML='.xml'



;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Reference Filter  & Reference Spectra              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--II.1-- Johnston V filter from Calar Alto  --------------------------------------------------
  refId  ='Generic/Johnson_UBVRIJHKL.V'
  fileRef='Generic.Johnson_UBVRIJHKL.V'
  if not file_test(root+fileRef+sufCSV) then begin
    refV = svoFilterSearch(ID=refId, dump=root+fileRef)
  endif else refV = svoFilterSearch(root+fileRef+sufXML)

 ;--II.2-- Reference Solar Spectrum  ----------------------------------------------------------
  dirSun = '/astrodata/Catalog/Taxonomy/SolarColors/templates/'
  readcol, dirSUN+'sun.csv', ws,fs,/silent
  sun={wave:ws*10, flux:fs,WavelengthUnit:'angstrom'}


 ;--II.3-- Reference Reflectivity of SSOs  ----------------------------------------------------
 ;--II.3.1-- Read Bus-DeMeo Asteroid Taxonomy  
  bdm = readBusDeMeo(/objects)
  bdmN = n_elements(bdm)

;--TAGGGG
;  bdmN = 5
  kbo={N:0}

 ;--II.3.2-- Read KBO Taxonomy from Merlin2017  -----------------------------------------------
;  dirKBO = '/astrodata/Catalog/Taxonomy/2016-Merlin/'
;  kboName = ['BB','BR','IR','RR']
;  nbKBO = n_elements(kboName)
;
;  readcol, dirKBO+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, dirKBO+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



;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---  III -- SSO Spectra in Physical Units                      -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--III.1-- SSO Structure  ---------------------------------------------------------------------
  eSSO = {name:'G2V', $
          wave:sun.wave, $
          flux:sun.flux, $
          dev:sun.flux*0.01, $
          WavelengthUnit:'angstrom'}

  nbSSO = 1 + bdmN
  sso = replicate(eSSO, nbSSO)

 ;--III.2-- Asteroid Spectra  ------------------------------------------------------------------
  for kS=0, bdmN-1 do begin
  
   ;--III.2.1-- Taxonomic Class
    sso[1+kS].name = bdm[kS].name

   ;--III.2.2-- Wavelength Range
    locWave=(*bdm[kS].wave)*1e4 
    minBDM = min(locWave, max=maxBDM)
    lowBDM  = where( sun.wave le minBDM, nbLow )
    highBDM = where( sun.wave ge maxBDM, nbHigh )
  
   ;--III.2.3-- Interpolation on Solar Spectrum
    newSpec = interpol( *bdm[kS].spec,                locWave, sun.wave )
    newLow  = interpol( *bdm[kS].spec + *bdm[kS].unc, locWave, sun.wave )
    newHigh = interpol( *bdm[kS].spec - *bdm[kS].unc, locWave, sun.wave )

   ;--III.2.4-- Clean Outside Taxonomy Definition
    newSpec[highBDM] = newSpec[highBDM[0]]
    newSpec[lowBDM]  = newSpec[lowBDM[nbLow-1]]
    newHigh[highBDM] = newHigh[highBDM[0]]
    newHigh[lowBDM]  = newHigh[lowBDM[nbLow-1]]
    newLow [highBDM] = newLow [highBDM[0]]
    newLow [lowBDM]  = newLow [lowBDM[nbLow-1]]

;    window, 2
;    cgPlot, locWave, bdm.class[kS].spec, thick=2, xRange=[0,5]*1e4, title=sso[1+kS].name
;    cgPlot, /over, sso[0].wave, newSpec, linestyle=2, color='red', thick=4
;    cgPlot, /over, sso[0].wave, newLow,  linestyle=2, color='Gray', thick=2
;    cgPlot, /over, sso[0].wave, newHigh, linestyle=2, color='Gray', thick=2

   ;--III.2.5-- Flux Density 
    sso[1+kS].flux = sun.flux * newSpec         
    sso[1+kS].dev  = sun.flux * (newHigh-newLow)

;    window, 2
;    cgPlot, sun.wave, sun.flux, thick=2, xRange=[0,2.5]*1e4, title=sso[kS+1].name
;    cgPlot, /over, sso[1+kS].wave, sso[1+kS].flux, linestyle=2, color='red', thick=2

   ;--III.2.6-- Write on Disk
    forprint, sun.wave, sso[1+kS].flux, abs(sso[1+kS].dev), $
              format='(E10.3,", ",E10.3,", ",E10.3)', /Silent, $
              textout=dirSSO+sso[1+kS].name+'.csv', $
              comment='Wavelength(nm),Flux,Deviation'
    

  endfor

  
 ;--III.3-- KBOs Spectra  ----------------------------------------------------------------------
;  for kS=0, kbo.N-1 do begin
;  
;   ;--III.3.1-- Taxonomic Class
;    sso[1+bdmN+kS].name = kbo.class[kS].name
;
;   ;--III.3.2-- Wavelength Range
;    locWave=kbo.wave*1e4 
;    minKBO = min(locWave, max=maxKBO)
;    lowKBO  = where( sun.wave le minKBO, nbLow )
;    highKBO = where( sun.wave ge maxKBO, nbHigh )
;  
;   ;--III.3.3-- Interpolation on Solar Spectrum
;    newSpec = interpol( kbo.class[kS].spec,                   locWave, sun.wave )
;    newLow  = interpol( kbo.class[kS].spec+kbo.class[kS].dev, locWave, sun.wave )
;    newHigh = interpol( kbo.class[kS].spec-kbo.class[kS].dev, locWave, sun.wave )
;
;   ;--III.3.4-- Clean Outside Taxonomy Definition
;    newSpec[highKBO] = newSpec[highKBO[0]]
;    newSpec[lowKBO]  = newSpec[lowKBO[nbLow-1]]
;    newHigh[highKBO] = newHigh[highKBO[0]]
;    newHigh[lowKBO]  = newHigh[lowKBO[nbLow-1]]
;    newLow [highKBO] = newLow [highKBO[0]]
;    newLow [lowKBO]  = newLow [lowKBO[nbLow-1]]
;
;;    window, 2
;;    cgPlot, locWave, kbo.class[kS].spec, thick=2, xRange=[0,5]*1e4, title=sso[1+bdmN+kS].name
;;    cgPlot, /over, sso[0].wave, newSpec, linestyle=2, color='red', thick=4
;;    cgPlot, /over, sso[0].wave, newLow,  linestyle=2, color='Gray', thick=2
;;    cgPlot, /over, sso[0].wave, newHigh, linestyle=2, color='Gray', thick=2
;
;   ;--III.3.5-- Flux Density 
;    sso[1+bdmN+kS].flux = sun.flux * newSpec         
;    sso[1+bdmN+kS].dev  = sun.flux * (newHigh-newLow)
;
;;    window, 2
;;    cgPlot, sun.wave, sun.flux, thick=2, xRange=[0,2.5]*1e4, title=sso[1+bdmN+kS].name
;;    cgPlot, /over, sso[1+bdmN+kS].wave, sso[1+bdmN+kS].flux, linestyle=2, color='red', thick=2
;
;   ;--III.3.6-- Write on Disk
;    forprint, sun.wave, sso[1+bdmN+kS].flux, abs(sso[1+bdmN+kS].dev), $
;              format='(E10.3,", ",E10.3,", ",E10.3)', /Silent, $
;              textout=dirSSO+sso[1+bdmN+kS].name+'.csv', $
;              comment='Wavelength(nm),Flux,Deviation'
;    
;
;  endfor


 ;--III.5-- Reference Magnitude in V  ----------------------------------------------------------
  magV = fltarr(nbSSO)
  headLine='Facility, Instrument, filterId, Band'
  for kSSO=0, nbSSO-1 do begin
    headLine+=', '+sso[kSSO].name
    magV[kSSO] = specToMag(sso[kSSO], refV, AB=AB, ST=ST, Vega=Vega)
  endfor

;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   IV -- Facility/Instrument                                -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  VminusF = fltarr(10,nbSSO)
  kConf=0
  for kSet=0, nbSet-1 do begin
    print, kSet+1,nbSet,set[kSet].facility, set[kSet].instrument, $
           format='(I2,"/",I2,2x,A-10,1x,A-10)'

   ;--IV.1-- Facility Directory  ---------------------------------------------------------------
   ;--IV.1.1-- Create Directory
    dirFAC = root+set[kSet].facility+'/'
    if not file_test(dirFAC,/Directory) then file_mkdir, dirFAC

   ;--IV.1.2-- Get Facility Information - Generic Case 
    if strCmp(set[kSet].facility, 'Generic') then begin
      facXML = dirFAC+set[kSet].instrument+sufXML
      facCSV = dirFAC+set[kSet].instrument+sufCSV
      if not file_test( facCSV, /Read ) then begin
        facInfo=svoFilterSearch(photSystem=set[kSet].instrument, dump=dirFAC+set[kSet].instrument)
      endif else facInfo=svoFilterSearch(facCSV)

   ;--IV.1.2-- Get Facility Information - Specific Facility
    endif else begin
      facXML = dirFAC+set[kSet].facility+sufXML
      facCSV = dirFAC+set[kSet].facility+sufCSV
      if not file_test( facCSV, /read ) then begin
        facInfo=svoFilterSearch(facility=set[kSet].facility, dump=dirFAC+set[kSet].facility)
      endif else facInfo=svoFilterSearch(facCSV)

    endelse

   ;--IV.1.3-- Validation of Results
    if size(facInfo,/Type) ne 8 then begin
      print, 'Facility cannot be found: '+set[kSet].facility
      stop
    endif


   ;--IV.2-- Instrument Directory  ------------------------------------------------------------
   ;--IV.2.1-- Select Generic from List
    if strCmp(set[kSet].facility, 'Generic') then begin
      pFilt = where( strCmp('Generic', facInfo.filterId, 7), nbFilt )
      if nbFilt eq 0 then begin 
        print, 'Generic '+set[kSet].instrument+' cannot be found.'
        stop
      endif
      

   ;--IV.2.1/B-- Select Instrument
    endif else begin

     ;--IV.2.1/B.1-- All Instruments
       if strCmp(set[kSet].instrument,'*') then begin
         nbFilt=n_elements(facInfo)
         pFilt=indgen(nbFilt)

     ;--IV.2.1/B.2-- Specific Instrument
      endif else begin
        pFilt = where( strCmp(set[kSet].instrument, facInfo.instrument), nbFilt )
        if nbFilt eq 0 then begin 
          print, 'Instrument '+set[kSet].instrument+' cannot be found in '+set[kSet].facility
          stop
        endif
      endelse
    endelse

   ;--IV.2.2-- Create Directory
    if strCmp(set[kSet].instrument,' ') then set[kSet].instrument=''
    dirINS = root+set[kSet].facility+'/'+set[kSet].instrument+'/'
    if not file_test(dirINS,/Directory) then file_mkdir, dirINS


   ;--IV.3-- Output Structures  ---------------------------------------------------------------
    colArr = fltarr(nbFilt, nbSSO)
    set[kSet].file = dirINS+ssoColor
    openW, idW, set[kSet].file, /Get_Lun
    printf, idW, headLine

   ;--IV.4-- Color Computation  ---------------------------------------------------------------
    for kF=0, nbFilt-1 do begin

     ;--IV.4.1-- Get Filter Curve from SVO Filter
      split=strSplit(facInfo[pFilt[kF]].filterId, '/', /Extract, count=nbSplit)
      if nbSplit eq 1 then filtName=facInfo[pFilt[kF]].filterId else filtName=split[1]
      filtXML = dirINS+filtName+sufXML
      filtCSV = dirINS+filtName+sufCSv
      if not file_test( filtXML, /Read ) then begin
        filtInfo=svoFilterSearch(ID=facInfo[pFilt[kF]].filterId, dump=dirINS+filtName)
      endif else filtInfo=svoFilterSearch(filtXML)


      instOut = set[kSet].instrument
      facOut  = set[kSet].facility
      bandOut = facInfo[pFilt[kF]].Band
      case set[kSet].facility of
         'SLOAN': instOut = 'SDSS'
         '2MASS': instOut = '2MASS'
         'Paranal': instOut = 'VISTA'
         'GAIA': begin
            instOut = 'Gaia'
            case facInfo[pFilt[kF]].filterId of
              'GAIA/GAIA2r.Gbp': bandOut='BP'
              'GAIA/GAIA2r.Grp': bandOut='RP'
              'GAIA/GAIA2r.G'  : bandOut='G'
            endcase
         end
         'Kepler': bandOut = 'K'
         else: begin
            instOut = set[kSet].instrument
            facOut  = set[kSet].facility
            bandOut = facInfo[pFilt[kF]].Band
         end
      endcase
      line= facOut+', '+instOut+', '+facInfo[pFilt[kF]].filterId+', '+bandOut


     ;--IV.4.2-- Compute Magnitude in Current Filter
      for kSSO=0, nbSSO-1 do begin
        magFilt = specToMag(sso[kSSO], filtInfo, AB=AB, ST=ST, Vega=Vega)
        colArr[kF,kSSO] = magV[kSSO] - magFilt
        line+= ', '+string(colArr[kF,kSSO],format='(F8.4)')

        VminusF[kConf,kSSO] = colArr[kF,kSSO]

      endfor
      printf, idW, line
      
      kConf++
    endfor    ;--IV.4-- Loop over Filters
    free_lun, idW

  endfor    ;--IV-- Loop over Sets



  openW, idW, root+ssoOut, /Get_Lun
  printf, idW, 'Number,Object,Class,u,g,r,i,z,Y,J,H,Ks'
  for kS=1, nbSSO-1 do begin
    line=string(bdm[kS-1].num, format='(I6)')+','+$
         string(bdm[kS-1].name,format='(A-12)')+','+$
         string(bdm[kS-1].class,format='(A-3)')
    sel=[0,1,2,3,4,6,7,8,9]
    nbSel=n_elements(sel)
    for kF=0, nbSel-1 do begin
      line += ','+string(VminusF[sel[kF],kS],format='(F8.4)')
    endfor
    printf, idW, line
  endfor
  free_lun, idW


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    V -- Concatenation of Colors                            -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  outFile = root+ssoRoot+'-'+nameSet+sufSys
  spawn, 'cp '+set[0].file+' '+outFile
  for kSet=1, nbSet-1 do spawn, "awk 'NR>1{print}' "+set[kSet].file+' >> '+outFile

end