; Takes Bus-DeMeo taxonomy average spectra and converts them to VISTA colors ; Returns flux and errors function taxo2vista, demeoTAXO, vistaDIR, plotTAXO=plotTAXO, dirWORK=dirWORK xRang = [0.3, 2.5] yRang = [0.5, 2.1] ;--------------------------------------------------------------------; ;--------------------------------------------------------------------; ;-- TAG -- 1 -- Initialization, Constants... ----; ;--------------------------------------------------------------------; ;--------------------------------------------------------------------; if not keyword_set(dirWORK) then dirWORK = vistaDIR if not keyword_set(filtwave) then filtwave = [1.020,0.877,1.252,1.645,2.147] nbClass=n_elements( demeoTAXO.class ) bdflux = fltarr(nbClass,5) ; hold mean flux for each color for each class bdfluxerr= fltarr(nbClass,5) ; hold mean error for each color for each class vistaTAXO = {Y: fltarr(nbClass), $ Z: fltarr(nbClass), $ J: fltarr(nbClass), $ H: fltarr(nbClass), $ K: fltarr(nbClass) } bdflux2 = fltarr(nbClass,5) ; hold mean flux for each color for each class bdfluxerr2= fltarr(nbClass,5) ; hold mean error for each color for each class vistaTAXO2 = vistaTAXO ;-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_Y.dat', wY, tY, /SILENT readcol, vistaDIR+'VISTA_Filters_at80K_forETC_Z.dat', wZ, tZ, /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 ;---- modify filter transmission to account for M1/M2/QE lonY = fltarr(10001) minY = min(wY, max=maxY) gY=where( wY ge minY and wY le 2500, nY) lonY(minY-1:minY+nY-2)=tY(gY) ; lonY(minY-1:maxY-1)=tY print, 'Y', minY, 2500 lonZ = fltarr(10001) minZ = min(wZ, max=maxZ) lonZ(minZ-1:maxZ-1)=tZ print, 'Z', minZ, maxZ lonJ = fltarr(10001) minJ = min(wJ, max=maxJ) gJ=where( wJ ge minJ and wJ le 2500, nJ) lonJ(minJ-1:minJ+nJ-2)=tJ(gJ) ; lonJ(minJ-1:maxJ-1)=tJ print, 'J', minJ, 2500 lonH = fltarr(10001) minH = min(wH, max=maxH) gH=where( wH ge minH and wH le 2500, nH) lonH(minH-1:minH+nH-2)=tH(gH) ; lonH(minH-1:maxH-1)=tH print, 'H', minH, 2500 lonK = fltarr(10001) minK = min(wK, max=maxK) gK=where( wK ge minK and wK le 2500, nK) lonK(minK-1:minK+nK-2)=tK(gK) ; lonK(minK-1:maxK-1)=tK print, 'K', minK, 2500 lonQE = fltarr(10001) minQE = min(wQE, max=maxQE) lonQE(minQE-1:maxQE-1)=tQE print, 'QE', minQE, maxQE ;-- ;--; lonM1 = fltarr(10001) ;--; minM1 = min(wM1, max=maxM1) ;--; gM1=where( wM1 ge minM1 and wM1 le 2500, nM1) ;--; lonM1(minM1-1:minM1+nM1-2)=tM1(gM1) ;--;; lonM1(minM1-1:maxM1-1)=tM1 ;--; print, 'M1', minM1, 2500 ;--; ;--; lonM2 = fltarr(10001) ;--; minM2 = min(wM2, max=maxM2) ;--; gM2=where( wM2 ge minM2 and wM2 le 2500, nM2) ;--; lonM2(minM2-1:minM2+nM2-2)=tM2(gM2) ;--;; lonM2(minM2-1:maxM2-1)=tM2 ;--; print, 'M2', minM2, 2500 ;--; ;-- ;-- lonYc = lonY*lonQE/100. lonZc = lonZ*lonQE/100. lonJc = lonJ*lonQE/100. lonHc = lonH*lonQE/100. lonKc = lonK*lonQE/100. lonW = indgen(10001) ;-- cgPlot, lonW, lonY, xr=[500,2600] ;-- cgPlot, /over, lonW, lonZ ;-- cgPlot, /over, lonW, lonJ ;-- cgPlot, /over, lonW, lonH ;-- cgPlot, /over, lonW, lonK ;-- cgplot, /over, lonW, lonQE, color='blue' ;--; cgplot, /over, lonW, lonM1, color='orange' ;--; cgplot, /over, lonW, lonM2, color='orange' ;-- cgPlot, /over, lonW, lonYc, linestyle=2 ;-- cgPlot, /over, lonW, lonZc, linestyle=2 ;-- cgPlot, /over, lonW, lonJc, linestyle=2 ;-- cgPlot, /over, lonW, lonHc, linestyle=2 ;-- cgPlot, /over, lonW, lonKc, linestyle=2 ;-- ;-- ;--; tZ2 = tZ * tQE ;--stop ;--; lonW /= 1.e3 ;-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 ) validZ2 = where( lonZc ne 0 ) validY2 = where( lonYc ne 0 ) validJ2 = where( lonJc ne 0 ) validH2 = where( lonHc ne 0 ) validK2 = where( lonKc ne 0 ) if keyword_set(plotTAXO) eq 1 then begin set_plot,'ps' device, xs=29.7, ys=21 device, filename=dirWORK+'vista-taxonomy.eps', /color, encapsulated=0 loadct, 0, /silent !X.margin=[4,0.5] !Y.margin=[2,2] multiplot, [5,5] !P.CHARSIZE = 0.8 !X.STYLE=1 & !Y.STYLE=1 !Y.TICKINTERVAL=.5 & !X.TICKINTERVAL=.2 !Y.TICKLEN=.03 !X.TICKLEN=.05 endif ;--------------------------------------------------------------------; ;--------------------------------------------------------------------; ;-- TAG -- 2 -- Compute Fluxes, Colors... ----; ;--------------------------------------------------------------------; ;--------------------------------------------------------------------; normaI = 2 cutTrans = 25. for kClass=0, nbClass-1 do begin specY = interpol( demeoTAXO.spec(kClass,*), demeoTAXO.wave, wY ) specZ = interpol( demeoTAXO.spec(kClass,*), demeoTAXO.wave, wZ ) specJ = interpol( demeoTAXO.spec(kClass,*), demeoTAXO.wave, wJ ) specH = interpol( demeoTAXO.spec(kClass,*), demeoTAXO.wave, wH ) specK = interpol( demeoTAXO.spec(kClass,*), demeoTAXO.wave, wK ) errY = interpol( demeoTAXO.err(kClass,*), demeoTAXO.wave, wY ) errZ = interpol( demeoTAXO.err(kClass,*), demeoTAXO.wave, wZ ) errJ = interpol( demeoTAXO.err(kClass,*), demeoTAXO.wave, wJ ) errH = interpol( demeoTAXO.err(kClass,*), demeoTAXO.wave, wH ) errK = interpol( demeoTAXO.err(kClass,*), demeoTAXO.wave, wK ) CASE normaI OF 0: BEGIN normaR = where( tY ge cutTrans ) normaV = mean(specY( normaR )) END 1: BEGIN normaR = where( tZ ge cutTrans ) normaV = mean(specZ( 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 saveS = demeoTAXO.spec(kClass,*) saveE = demeoTAXO.err(kClass,*) demeoTAXO.spec(kClass,*) /= normaV demeoTAXO.err(kClass,*) /= normaV specY /= normaV specZ /= normaV specJ /= normaV specH /= normaV specK /= normaV ;--2.1-- Convert DeMeo Spectra into VISTA Colors (Normalized to R) fluxes=[ total( specY(validY)*tY(validY) ) / total(tY(validY)), $ total( specZ(validZ)*tZ(validZ) ) / total(tZ(validZ)), $ 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.Y(kClass) = colors(0) vistaTAXO.Z(kClass) = colors(1) vistaTAXO.J(kClass) = colors(2) vistaTAXO.H(kClass) = colors(3) vistaTAXO.K(kClass) = colors(4) bdflux(kClass,*) = [ 10.^( -0.4* (vistaTAXO.Y(kClass)) ), $ 10.^( -0.4* (vistaTAXO.Z(kClass)) ), $ 10.^( -0.4* (vistaTAXO.J(kClass)) ), $ 10.^( -0.4* (vistaTAXO.H(kClass)) ), $ 10.^( -0.4* (vistaTAXO.K(kClass)) ) ] ;--2.2-- Associate DeMeo Deviation into VISTA Colors bdfluxerr(kClass,*)=[ mean(errY(validY)), $ mean(errZ(validZ)), $ mean(errJ(validJ)), $ mean(errH(validH)), $ mean(errK(validK)) ] ;------------------------------test specY2 = interpol( saveS, demeoTAXO.wave, lonW ) specZ2 = interpol( saveS, demeoTAXO.wave, lonW ) specJ2 = interpol( saveS, demeoTAXO.wave, lonW ) specH2 = interpol( saveS, demeoTAXO.wave, lonW ) specK2 = interpol( saveS, demeoTAXO.wave, lonW ) errY2 = interpol( saveE, demeoTAXO.wave, lonW ) errZ2 = interpol( saveE, demeoTAXO.wave, lonW ) errJ2 = interpol( saveE, demeoTAXO.wave, lonW ) errH2 = interpol( saveE, demeoTAXO.wave, lonW ) errK2 = interpol( saveE, demeoTAXO.wave, lonW ) normaR2 = where( lonJc ge cutTrans ) normaV2 = mean(specJ2( normaR2 )) specY2 /= normaV2 specZ2 /= normaV2 specJ2 /= normaV2 specH2 /= normaV2 specK2 /= normaV2 ;; cgPlot, demeoTAXO.wave, saveS, xr=[0.5,2.5], thick=2 ;; cgPlot, /over, lonW, specY2, color='cornflower blue', thick=2 ;--2.1-- Convert DeMeo Spectra into VISTA Colors (Normalized to R) fluxes2=[ total( specY2(validY2)*lonYc(validY2) ) / total(lonYc(validY2)), $ total( specZ2(validZ2)*lonZc(validZ2) ) / total(lonZc(validZ2)), $ total( specJ2(validJ2)*lonJc(validJ2) ) / total(lonJc(validJ2)), $ total( specH2(validH2)*lonHc(validH2) ) / total(lonHc(validH2)), $ total( specK2(validK2)*lonKc(validK2) ) / total(lonKc(validK2)) ] ;; cgplot, /over, lonW, specY2*lonYc/100 ;; cgplot, /over, lonW, specZ2*lonZc/100 ;; cgplot, /over, lonW, specJ2*lonJc/100 ;; cgplot, /over, lonW, specH2*lonHc/100 ;; cgplot, /over, lonW, specK2*lonKc/100 colors2 = -2.5*alog10( fluxes2/fluxes2(normaI) ) vistaTAXO2.Y(kClass) = colors2(0) vistaTAXO2.Z(kClass) = colors2(1) vistaTAXO2.J(kClass) = colors2(2) vistaTAXO2.H(kClass) = colors2(3) vistaTAXO2.K(kClass) = colors2(4) ;; cgPlot, /over, filtwave, fluxes2, psym='Filled circle' bdflux2(kClass,*) = [ 10.^( -0.4* (vistaTAXO2.Y(kClass)) ), $ 10.^( -0.4* (vistaTAXO2.Z(kClass)) ), $ 10.^( -0.4* (vistaTAXO2.J(kClass)) ), $ 10.^( -0.4* (vistaTAXO2.H(kClass)) ), $ 10.^( -0.4* (vistaTAXO2.K(kClass)) ) ] ;--2.2-- Associate DeMeo Deviation into VISTA Colors bdfluxerr2(kClass,*)=[ mean(errY2(validY2)), $ mean(errZ2(validZ2)), $ mean(errJ2(validJ2)), $ mean(errH2(validH2)), $ mean(errK2(validK2)) ] ;--- show difference between colours with and witout QE ;print, reform(bdflux(kClass,*)/bdflux2(kClass,*)) ;print, reform(colors-colors2), format='(5(F7.4,2x))' if keyword_set(plotTAXO) eq 1 then begin ;--------------------------------------------------------------------------- ;--2.3-- TaxoPlot ;--2.3.1-- Preparation plot, demeoTAXO.wave, demeoTAXO.spec(kClass,*), $ xr=xRang, yr=yRang, xst=1,yst=1,/NODATA, $ xtickint=0.5, xminor=5 ;--2.3.2-- Taxonomic class range range=where( demeoTAXO.wave ge 0.35 and demeoTAXO.wave le 2.4 ) xClass = [ reform(demeoTAXO.wave(range)), $ reform(reverse(demeoTAXO.wave(range)))] yClass = [ reform(demeoTAXO.spec(kClass,range)-demeoTAXO.err(kClass,range)), $ reform(demeoTAXO.spec(kClass,reverse(range))+demeoTAXO.err(kClass,reverse(range))) ] polyfill, xClass, yClass, /DATA, color=200 oplot, demeoTAXO.wave, demeoTAXO.spec(kClass,*) oplot, xClass, yClass, color=130 ;--2.3.3-- VISTA Color for Taxonomy Classes define_psym, 1, symbol=11 oploterror, filtwave, bdflux(kClass,*), bdfluxerr(kClass,*), psym=3 oplot, filtwave, bdflux(kClass,*), psym=8 ;; cgPlot, /over, filtwave, bdflux2(kClass,*), psym='Filled square', color='red' ; OPLOT, wY, specY+0.4, psym=3, color=80 , thick=2 ; OPLOT, wZ, specZ+0.2, psym=3, color=230, thick=2 ; OPLOT, wJ, specJ+0.0, psym=3, color=130, thick=2 ; OPLOT, wH, specH-0.2, psym=3, color=230, thick=2 ; OPLOT, wK, specK-0.4, psym=3, color=80 , thick=2 ;--2.3.4-- Class Name xyouts,0.5,1.5, demeoTAXO.class(kclass), charsize=1.5, charthick=2 ; oplot, wY, tY/50. ; oplot, wZ, tZ/50. ; oplot, wJ, tJ/50. ; oplot, wH, tH/50. ; oplot, wK, tK/50. multiplot endif endfor if keyword_set(plotTAXO) eq 1 then begin plot, wY, tY, /DEVICE, $ xr=xRang, yr=[0,100], xst=1,yst=1, $ title='Filter Response Curve', $ position=[24000,700,29600,4000], $ ytickint=10 oplot, wZ, tZ oplot, wJ, tJ oplot, wH, tH oplot, wK, tK cgPlot, /over, wQE, tQE, color='Gray' ; oplot, wM1, tM1 ; oplot, wM2, tM2 ; xyouts, 0.21, 0.5, '0.5' ; we are not printing taxo-vista.txt unless we plot fmt='(10(F8.4,3x))' filen=dirWORK+'taxo-vista.txt' forprint,bdflux(*,0),bdfluxerr(*,0),$ bdflux(*,1),bdfluxerr(*,1),$ bdflux(*,2),bdfluxerr(*,2),$ bdflux(*,3),bdfluxerr(*,3),$ bdflux(*,4),bdfluxerr(*,4), format=fmt, textout=filen,/nocomment,/silent multiplot,/reset xyouts, 15000, 20500, 'Bus-DeMeo Taxonomy Converted into VISTA Colors', /DEVICE, align=0.5, chars=1.5, chart=2 device,/close set_plot,'X' endif bd = {ref: bdflux, $ err: bdfluxerr } return, bd end demeo = taxoRead('/obs/bcarry/data/mining/neas/bus-demeo/meanspectra.tab' ) visDIR = '/obs/bcarry/data/mining/neas/VISTA/' filtwave = [1.020,0.877,1.252 ,1.645 ,2.147] taxo = taxo2vista( demeo, visDIR, filtWave,/plot) ; taxo = taxo2vista( demeo, visDIR, filtWave ) end