;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- General Configuration ---------------------------------------------------------------; ; cwd = '/home/bcarry/work/data/lcth/' cwd = '/home/bcarry/data/lc/lcth/' debug=0 ;--II.2-- Solar-like Spectra Parameters ------------------------------------------------------; sunTemp = 5778. fileG2V = '/astrodata/Catalog/Taxonomy/SolarColors/templates/lte058-4.5-0.0a+0.0.BT-NextGen.7.dat.csv' fileSun = '/astrodata//Catalog/Taxonomy/SolarColors/templates/sun.csv' ;--II.3-- Wavelength Limits ------------------------------------------------------------------; vBand = [0.463, 0.639] ;-micron limitBlue = 0.450 ;-micron limitRed = 2.450 ;-micron rangeV = [0.463, 0.639] ;-micron rangeTaxo = [0.450, 2.450] ;-micron rangeWave = [ 1e-5, 1e2 ] ;-micron limitFrac = 0.01 ;--II.4-- Read Bus-DeMeo Taxonomy ------------------------------------------------------------; bdm = readBusDeMeo() ;--II.X-- xxxxxxx ----------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Solar-like Spectrum -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Blackbody Function at 5778 K -------------------------------------------------------; ;--II.1.1-- Create a Blackbody with Many Points at Max Energy nbWave = 1000 waveB = rangeWave[0] + (rangeTaxo[0]-rangeWave[0]) * findgen(nbWave)/(nbWave-1) nbWave = 10000 waveV = rangeTaxo[0] + (rangeTaxo[1]-rangeTaxo[0]) * findgen(nbWave)/(nbWave-1) nbWave = 1000 waveR = rangeTaxo[1] + (rangeWave[1]-rangeTaxo[1]) * findgen(nbWave)/(nbWave-1) wave= [waveB,waveV,waveR] BB = {wave: wave, flux:planck( wave*1e4, sunTemp)} ;-ergs/cm2/s/A ; cgPlot, BB.wave, BB.flux, /ylog, /xlog, xRange=[0.1, 5], yRange=[1e5,1e8] ;--II.1.2-- Interpolate Bus-DeMeo Spectra bdmBB = fltarr(n_elements(BB.wave),bdm.N) for kC=0, bdm.N-1 do begin iBDM = where( BB.wave ge rangeTaxo[0] and BB.wave le rangeTaxo[1], nbIndex ) interp = interpol( bdm.class[kC].spec, bdm.wave, BB.wave[iBDM] ) bdmBB[iBDM,kC] = interp low = where( BB.wave lt rangeTaxo[0] ) bdmBB[low,kC] = bdmBB[iBDM[0],kC] high = where( BB.wave gt rangeTaxo[1] ) bdmBB[high,kC] = bdmBB[iBDM[-1],kC] bdmBB[*,kC] *= BB.flux ; cgPlot, BB.wave, bdmBB[*,kC], /overplot, color=bdm.class[kC].color endfor ;--II.2-- Read G2V star Synthetic Spectrum ---------------------------------------------------; ;--II.2.1-- G2V Spectrum from File readcol, fileG2V, w, f, delimiter=',', /Silent g2v = {wave: w*1e-4, flux:f } ; cgPlot, g2v.wave, g2v.flux, /ylog, /xlog, xRange=[0.1, 5], yRange=[1e5,1e8] ;--II.1.2-- Interpolate Bus-DeMeo Spectra bdmG2 = fltarr(n_elements(g2v.wave),bdm.N) for kC=0, bdm.N-1 do begin iBDM = where( g2v.wave ge rangeTaxo[0] and g2v.wave le rangeTaxo[1], nbIndex ) interp = interpol( bdm.class[kC].spec, bdm.wave, g2v.wave[iBDM] ) bdmG2[iBDM,kC] = interp low = where( g2v.wave lt rangeTaxo[0] ) bdmG2[low,kC] = bdmG2[iBDM[0],kC] high = where( g2v.wave gt rangeTaxo[1] ) bdmG2[high,kC] = bdmG2[iBDM[-1],kC] bdmG2[*,kC] *= g2v.flux ; cgPlot, g2v.wave, bdmG2[*,kC], /overplot, color=bdm.class[kC].color endfor ;--II.3-- Real Solar Spectrum ----------------------------------------------------------------; ;--II.3.1-- Solar Spectrum from File readcol, fileSun, w, f, delimiter=',', /Silent sun = {wave: w*1e-3, flux:f } ; cgPlot, sun.wave, sun.flux, /ylog, /xlog, xRange=[0.1, 5], yRange=[0.01,2.5] ;--II.1.2-- Interpolate Bus-DeMeo Spectra bdmSun = fltarr(n_elements(sun.wave),bdm.N) for kC=0, bdm.N-1 do begin iBDM = where( sun.wave ge rangeTaxo[0] and sun.wave le rangeTaxo[1], nbIndex ) interp = interpol( bdm.class[kC].spec, bdm.wave, sun.wave[iBDM] ) bdmSun[iBDM,kC] = interp low = where( sun.wave lt rangeTaxo[0] ) bdmSun[low,kC] = bdmSun[iBDM[0],kC] high = where( sun.wave gt rangeTaxo[1] ) bdmSun[high,kC] = bdmSun[iBDM[-1],kC] bdmSun[*,kC] *= sun.flux ; cgPlot, sun.wave, bdmSun[*,kC], /overplot, color=bdm.class[kC].color endfor ;--II.X-- Quality Control Figure -------------------------------------------------------------; ; cgPS_open, Filename=cwd+'solarSpectrum.eps', /Metric, /Decomposed, /Portrait, $ ; xSize=30, ySize=20, language_level=2, /quiet ; ; cgPlot, g2v.wave, g2v.flux, /xLog, /yLog, $ ; xTitle='Wavelength ('+cgGreek('mu')+'m)', $ ; yTitle='Flux (ergs/cm!U2!N/s/A)', $ ; xRange=[0.01,1e3], $ ; yRange=[1e-6,1e8] ; cgPlot, /OverPlot, g2v.wave, g2v.flux, color='red' ; cgPlot, /OverPlot, BB.wave, BB.flux ; ; cgPS_close, /png, Delete_PS=0 ;stop ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Compute Energy Bins -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Setup Models & Bands --------------------------------------------------------------; nbModel = 3 nbBand = 4 energy = fltarr( nbModel, bdm.N+1, nbBand ) emptySED = {bot:0., top:0., frac:0., N:0L, taxo:fltarr(bdm.N) } SED = replicate(emptySED, nbModel, 10000L) for kMod=0, nbModel-1 do begin ;--III.2-- Details for Models --------------------------------------------------------------; case kMod of ;--III.2.1-- Simple Blackbody 0: begin print, 'Simple Blackbody' fileOut='BB.csv' locWave = BB.wave locFlux = BB.flux locBDM = bdmBB end ;--III.2.2-- G2V Synthetic spectrum 1: begin print, 'G2V Synthetic spectrum' fileOut='G2V.csv' locWave = g2v.wave locFlux = g2v.flux locBDM = bdmG2 end ;--III.2.3-- Solar Spectrum 2: begin print, 'Solar spectrum' fileOut='sun.csv' locWave = sun.wave locFlux = sun.flux locBDM = bdmSun end endcase totFlux = total(locFlux) totBDM = total(locBDM,1) ;--III.3-- Summary Information for Large Spectral Bands ----------------------------------; for kBand=0, nbBand-1 do begin ;--III.3.1-- Details for Bands ---------------------------------------------------------; case kBand of ;--III.3.1/A-- From 0 to Blue 0: begin locLabel= 'Blue-end ' locMin = rangeWave[0] locMax = rangeTaxo[0] end ;--III.3.2/B-- DeMeo Wavelength Range 1: begin locLabel= 'DeMeo range ' locMin = rangeTaxo[0] locMax = rangeTaxo[1] end ;--III.3.2/C-- Johnston V Filter 2: begin locLabel= 'Johnston V ' locMin = rangeV[0] locMax = rangeV[1] end ;--III.3.3/D-- From Red to infinity 3: begin locLabel= 'Red-end ' locMin = rangeTaxo[1] locMax = rangeWave[1] end endcase ;--III.3.2-- Compute Energy Fraction in the Band ---------------------------------------; inBand = where( locWave ge locMin and locWave lt locMax ) energy[kMod,0,kBand] = total( locFlux[inBand] ) / totFlux for kC=0, bdm.N-1 do energy[kMod,kC+1,kBand] = total( locBDM[inBand,kC] ) / totBDM[kC] print, ' '+locLabel+ string( 100*energy[kMod,0,kBand],format='(F4.1)')+' %' endfor print, ' Total '+ string( 100*total(energy[kMod,0,[0,1,3]]),format='(F6.1)')+' %' ; for kC=0,bdm.N-1 do $ ; print, ' '+string(bdm.class[kC].name,format='(A3)')+'-type : '+string( 100*total(energy[kMod,kC+1,[0,1,3]]),format='(F6.1)')+' %' ;--- bdm gives 100% too ; cgPlot, locWave, locFlux, /ylog, /xlog, xRange=[0.01, 1e3], yRange=[1e-4,1e8] ;--III.4-- Detailed SED ------------------------------------------------------------------; kBin=0L ;--III.4.1-- From Minimum Wavelength to Bus-DeMeo Interval -------------------------------; ;--III.4.1/A-- Wavelength Range & Index locMin = rangeWave[0] locMax = rangeTaxo[0] index = where( locWave ge locMin and locWave lt locMax, nbEval ) ;--III.4.1/B-- Bookkeep Values SED[kMod,kBin].N = nbEval SED[kMod,kBin].bot = locMin SED[kMod,kBin].top = locMax SED[kMod,kBin].frac = energy[kMod,0,0] for kC=0,bdm.N-1 do SED[kMod,kBin].taxo[kC] = energy[kMod,kC+1,0] kBin++ ; cgPlot, /OverPlot, locWave[index], locFlux[index], color='Navy', psym='Filled Circle' print, '4.1: '+string(nbEval) ;--III.4.2-- From Minimum Bus-DeMeo to Johnston V ----------------------------------------; ;--III.4.2/A-- Wavelength Range & Index locMin = rangeTaxo[0] locMax = rangeV[0] index = where( locWave ge locMin and locWave lt locMax, nbEval ) ; cgPlot, /OverPlot, locWave[index], locFlux[index], color='Cornflower Blue', psym='Filled Circle' print, '4.2: '+string(nbEval) ;--III.4.2/B-- Loop to get 1% Total Energy Bins iMin = min( index, max=iMax ) kL=iMin kH=iMin+1 while kH le iMax do begin tVal = total( locFlux[kL:kH] ) / totFlux if tVal ge limitFrac then begin ;-- Got Above 1% - Make One Step Back kH-- tVal = total( locFlux[kL:kH] ) / totFlux ;-- Bookkeep Values SED[kMod,kBin].N = kH-kL SED[kMod,kBin].bot = locWave[kL] SED[kMod,kBin].top = locWave[kH] SED[kMod,kBin].frac = tVal for kC=0,bdm.N-1 do SED[kMod,kBin].taxo[kC] = total( locBDM[kL:kH,kC] ) / totBDM[kC] ;-- Setup New Step kL=kH kBin++ endif kH++ endwhile ;--III.4.3-- Inside Johnston V Filter ----------------------------------------------------; ;--III.4.3/A-- Wavelength Range & Index locMin = rangeV[0] locMax = rangeV[1] index = where( locWave ge locMin and locWave lt locMax, nbEval ) ;--III.4.3/B-- Bookkeep Values SED[kMod,kBin].N = nbEval SED[kMod,kBin].bot = locMin SED[kMod,kBin].top = locMax SED[kMod,kBin].frac = energy[kMod,0,2] for kC=0,bdm.N-1 do SED[kMod,kBin].taxo[kC] = energy[kMod,kC+1,2] kBin++ ; cgPlot, /OverPlot, locWave[index], locFlux[index], color='Green', psym='Filled Circle' print, '4.3: '+string(nbEval) ;--III.4.4-- From Johnston V to Maximum Bus-DeMeo ----------------------------------------; ;--III.4.4/A-- Wavelength Range & Index locMin = rangeV[1] locMax = rangeTaxo[1] index = where( locWave ge locMin and locWave lt locMax, nbEval ) ; cgPlot, /OverPlot, locWave[index], locFlux[index], color='Orange', psym='Filled Circle' print, '4.4: '+string(nbEval) ;--III.4.4/B-- Loop to get 1% Total Energy Bins iMin = min( index, max=iMax ) kL=iMin kH=iMin+1 while kH le iMax do begin tVal = total( locFlux[kL:kH] ) / totFlux if tVal ge limitFrac then begin ;-- Got Above 1% - Make One Step Back kH-- tVal = total( locFlux[kL:kH] ) / totFlux ;-- Bookkeep Values SED[kMod,kBin].N = kH-kL SED[kMod,kBin].bot = locWave[kL] SED[kMod,kBin].top = locWave[kH] SED[kMod,kBin].frac = tVal for kC=0,bdm.N-1 do SED[kMod,kBin].taxo[kC] = total( locBDM[kL:kH,kC] ) / totBDM[kC] ;-- Setup New Step kH++ kL=kH kBin++ endif kH++ endwhile ;--III.4.5-- From Maximum Bus-DeMeo to Maximum Wavelength --------------------------------; ;--III.4.5/A-- Wavelength Range & Index locMin = rangeTaxo[1] locMax = rangeWave[1] index = where( locWave ge locMin and locWave lt locMax, nbEval ) ;--III.4.5/B-- Bookkeep Values SED[kMod,kBin].N = nbEval SED[kMod,kBin].bot = locMin SED[kMod,kBin].top = locMax SED[kMod,kBin].frac = energy[kMod,0,3] for kC=0,bdm.N-1 do SED[kMod,kBin].taxo[kC] = energy[kMod,kC+1,3] kBin++ ; cgPlot, /OverPlot, locWave[index], locFlux[index], color='Red', psym='Filled Circle' print, '4.5: '+string(nbEval) ;--III.4.6-- Store Results on Hard Drive -------------------------------------------------; ;--III.4.6.1-- Header and Open File header='Lower wavelength, Upper wavelength, Evaluation points, Fraction of energy' for kC=0, bdm.N-1 do header+=', '+bdm.class[kC].name openW, idEx, cwd+fileOut, /Get_Lun printf, idEx, header sep=', ' for kP=0, kBin-1 do begin ;--III.4.6.2-- Wavelength Range line = string(SED[kMod,kP].bot, format='(F10.6)')+sep+$ string(SED[kMod,kP].top, format='(F10.6)')+sep ;--III.4.6.3-- Evaluation points and Fraction of Energy line += string(SED[kMod,kP].N, format='(I6)')+sep+$ string(SED[kMod,kP].frac*100, format='(F10.6)') ;--III.4.6.4-- Fraction of Energy for Each Taxonomic Class for kC=0, bdm.N-1 do line += sep+string(SED[kMod,kP].taxo[kC]*100, format='(F10.6)') ;--III.4.6.5-- Write the line printf, idEx, line endfor free_lun, idEX ; forprint, SED[kMod,*].bot, SED[kMod,*].top, SED[kMod,*].N, SED[kMod,*].frac*100, $ ; SED[kMod,*].taxo[0]*100, $ ; textout=cwd+fileOut, /Silent, $ ; comment='Lower wavelength, Upper wavelength, Evaluation points, Fraction of energy', $ ; format='(F10.6,", ",F10.6,", ",I6,", ",F10.6,", ",F10.6)', $ ; subset=indgen(kBin) print, ' Total energy (%): '+strTrim(string(total(SED[kMod,0:kBin-1].frac)*100,format='(F7.4)'),2) for kC=0, bdm.N-1 do print, ' '+string(bdm.class[kC].name,format='(A3)')+'-type (%): '+$ strTrim(string(total(SED[kMod,0:kBin-1].taxo[kC])*100,format='(F9.4)'),2) print, '' print, '' endfor stop sunFlux = total(locFlux) print, sunFlux ; wMu = wave*1e-4 wMu = bb.wave ; cgWindow ; cgPlot, wMu, BB, /xLog, /yLog, /AddCmd, $ ; xtitle='Wavelength (micron)', ytitle='Flux (ergs/cm2/s/A)', $ ; xr=[0.01,10.], yr=[1e4,1e8] ;----Far red wMin=limitRed ;-micron wMax=10.000 ;-micron wAmin = wMin ; *1e4 wAmax = wMax ; *1e4 inFR = where( bb.wave ge wAmin and bb.wave le wAmax, nbFR) BB_FR = total( BB.flux[inFR] ) print, 'FarRed (%):', 100.*BB_FR/sunFlux, nbFR farRed= 100.*BB_FR/sunFlux ; cgPlot, /over, wMu(inFR), BB(inFR), color='Black', thick=2, /AddCmd ;----DeMeo band wMin=limitBlue;-micron wMax=limitRed ;-micron wAmin = wMin ; *1e4 wAmax = wMax ; *1e4 inBD = where( bb.wave ge wAmin and bb.wave le wAmax, nbBD) BB_BD = total( BB.flux[inBD] ) print, 'DeMeo (%):', 100.*BB_BD/sunFlux, nbBD demeo=100.*BB_BD/sunFlux ; cgPlot, /over, wMu(inBD), BB(inBD), color='cornflower blue', thick=2, /AddCmd ;----Blue wMin=0.01 ;-micron wMax=limitBlue;-micron wAmin = wMin ; *1e4 wAmax = wMax ; *1e4 inBlue = where( bb.wave ge wAmin and bb.wave le wAmax, nbBlue) BB_Blue = total(BB.flux[inBlue]) print, 'Blue (%):', 100.*BB_Blue/sunFlux, nbBlue blue = 100.*BB_Blue/sunFlux ; cgPlot, /over, wMu(inBlue), BB(inBlue), color='Green', thick=4, /AddCmd print, 'Total (%):', blue+demeo+farRed, (1-(blue+demeo+farRed)/100)*1e6, ' ppm' print, '' ;----V band wMin=vBand[0] ;-micron wMax=vBand[1] ;-micron wAmin = wMin ; *1e4 wAmax = wMax ; *1e4 inV = where( bb.wave ge wAmin and bb.wave le wAmax, nbV) BB_V = total(BB.flux[inV]) print, 'V Band (%):', 100.*BB_V/sunFlux, nbV V = 100.*BB_V/sunFlux ; cgPlot, /over, wMu(inV), BB(inV), color='red', thick=4, /AddCmd ;----miss band ; wMin=0.200 ;-micron ; wMax=0.400 ;-micron ; ; ; wAmin = wMin*1e4 ; wAmax = wMax*1e4 ; inMiss = where( wave ge wAmin and wave le wAmax, nbMiss) ; ; BB_Miss = total(BB(inMiss)) ; ; print, '0.2-0.4 (%):', 100.*BB_Miss/sunFlux, nbMiss ; miss = 100.*BB_Miss/sunFlux ; ; cgPlot, /over, wMu(inMiss), BB(inMiss), color='Green', thick=4, /AddCmd ; cgDelete print, '' stop ;------energy bins ;----in the red eArr = 0 lArr = vBand[1] vArr = 100.*BB_V/sunFlux intArr = vBand lLeft = max(inV)+1 lTest = lLeft+1 lEnd = max(inBD) while lTest le lEnd do begin testE = total(BB.flux[lLeft:lTest]) ratio = testE / sunFlux ; print, ratio if ratio ge 0.01 then begin eArr = [eArr, testE] lArr = [lArr, wMu[lTest]] vArr = [vArr, ratio*100.] intArr=[ [intArr], [wMu[lLeft],wMu[lTest]] ] lLeft=lTest+1 lTest++ endif lTest++ endwhile vArr = [vArr, farRed] lArr = [lArr, 500.] intArr=[ [intArr], [limitRed,500] ] ; forprint, intArr[0,*], intArr[1,*], vArr, textout=2 ;----in the blue eArr2 = 0 lArr2 = vBand[0] vArr2 = 100.*BB_V/sunFlux intArr2 = vBand lRight = min(inV)-1 lTest = lRight-1 lEnd = min(inBD) while lTest ge lEnd do begin testE = total(BB.flux[lTest:lRight]) ratio = testE / sunFlux if ratio ge 0.01 then begin eArr2 = [eArr2, testE] lArr2 = [lArr2, wMu[lTest]] vArr2 = [vArr2, ratio*100.] intArr2=[ [intArr2], [wMu[lTest],wMu[lRight]] ] lRight=lTest-1 lTest-- endif lTest-- endwhile nbBand = n_elements(lArr2) if nbBand ge 3 then begin intArr2=reverse( intArr2[*,1:nbBand-1], 2) eArr2 = reverse( eArr2[1:nbBand-1] ) lArr2 = reverse( lArr2[1:nbBand-1] ) vArr2 = reverse( vArr2[1:nbBand-1] ) endif else begin intArr2=intArr2[*,1:nbBand-1] eArr2 = eArr2[1:nbBand-1] lArr2 = lArr2[1:nbBand-1] vArr2 = vArr2[1:nbBand-1] endelse ; forprint, intArr2[0,*], intArr2[1,*], vArr2, textout=2 vArr2 = [blue, vArr2] lArr2 = [0.01,limitBlue] intArr2=[ [0.01,limitBlue], [intArr2] ] lArr = [lArr2,lArr] eArr = [eArr2,eArr] vArr = [vArr2,vArr] intArr=[ [intArr2], [intArr] ] nbBand = n_elements(lArr) forprint, intArr[0,*], intArr[1,*], vArr, textout=2, format='(F9.5," , ",F9.5," , ",F9.5)' print, nbBand stop stop ;top ;------asteroid reflectance ; class = ['S' ,'C' ] ; pv = [0.20,0.05] ; color = ['red','gray'] ; ; nbClass=2 root='/home/bcarry/work/data/lcth/' bd='busdemeo-meanspectra.csv' bd = readBusDeMeo( ) ; bd = readBusDeMeo( root+bd ) nbClass=n_elements(bd.class) ; range = where( wMu ge 0.4 and wMu le 2.5 ) ; NN=100 ; wave = 0.425 + (2.475-0.425)*findgen(NN)/(NN-1) refl = fltarr(nbClass, nbBand) midW = fltarr(nbBand) ; cgWindow ; ; cgPlot, lArr,/AddCmd ; cgPlot, lArr,/AddCmd, /over, psym='filled circle' ; ; cgPlot, bd.wave, /over,/AddCmd, color='grey' ; cgPlot, bd.wave,/AddCmd, /over, psym='filled circle', color='grey' ; cgDelete track=intarr(nbClass,nbBand) for kB=1, nbBand-2 do begin ; print, '' ; print, kB, intArr[*,kB], format='(2x,I2,4x,F5.3,3x,F5.3)' left = where( bd.wave le intArr[0,kB] ) lL = max( left ) left = where( bd.wave ge intArr[0,kB] ) lU = min( left ) right = where( bd.wave le intArr[1,kB] ) rL = max( right ) right = where( bd.wave ge intArr[1,kB] ) rU = min( right ) ; print, bd.wave[lL], bd.wave[lU], format='(4x,4x,F5.3,3x,F5.3)' ; print, bd.wave[rL], bd.wave[rU], format='(4x,4x,F5.3,3x,F5.3)' midW[kB] = mean( intArr[*,kB] ) for kC=0, nbClass-1 do begin if (lL eq rL) and (lU eq rU) then begin track[kC,kB]=1 ind = [lL,lU] locW = bd.wave[ ind ] locS = reform(bd.class[kC].spec[ind]) ; locS = reform(bd.spec[kC, ind ]) if locS[0] eq locS[1] then begin a=0 b=locS[0] endif else a = regress( locW, locS, const=b) refl[kC,kB] = a*midW[kB]+b if debug eq 1 then begin xx = [Ll,lU,rL,rU] cgPlot, bd.wave[xx], reform(bd.class[kC].spec[ xx ]) cgPlot, bd.wave[xx], reform(bd.class[kC].spec[ xx ]),/over, psym='filled square' ; cgPlot, bd.wave[xx], reform(bd.spec[kC, xx ]) ; cgPlot, bd.wave[xx], reform(bd.spec[kC, xx ]),/over, psym='filled square' cgPlot, intarr[*,kB],[0,0], thick=10,/over cgPlot, /over, midW[kB], refl[kC,kB], psym='filled circle', symsize=3 endif endif else begin track[kC,kB]=2 if debug eq 1 then begin xx = [Ll,lU,rL,rU] cgPlot, bd.wave[xx], reform(bd.class[kC].spec[ xx ]) cgPlot, bd.wave[xx], reform(bd.class[kC].spec[ xx ]),/over, psym='filled square' ; cgPlot, bd.wave[xx], reform(bd.spec[kC, xx ]) ; cgPlot, bd.wave[xx], reform(bd.spec[kC, xx ]),/over, psym='filled square' cgPlot, intarr[*,kB],[0,0], thick=10,/over endif ;-fraction of bin -- Left ind = [lL,lU] locW = bd.wave[ ind ] locS = reform(bd.class[kC].spec[ind]) ; locS = reform(bd.spec[kC, ind ]) if locS[0] eq locS[1] then begin a=0 b=locS[0] endif else a = regress( locW, locS, const=b) xArea = [ intArr[0,kB], bd.wave[[lU,lU]] , intArr[0,kB] ] yArea = [ 0,0, bd.class[kC].spec[lU], a*intArr[0,kB]+b ] ; yArea = [ 0,0, bd.spec[kC,lU], a*intArr[0,kB]+b ] aL = poly_area( xArea,yArea) sL = bd.wave[lU]-intArr[0,kB] if debug eq 1 then cgPlot, /over, xArea, yArea, thick=2 ;-fraction of bin -- Right ind = [rL,rU] locW = bd.wave[ ind ] locS = reform(bd.class[kC].spec[ind]) ; locS = reform(bd.spec[kC, ind ]) if locS[0] eq locS[1] then begin a=0 b=locS[0] endif else a = regress( locW, locS, const=b) xArea = [bd.wave[rL], intArr[1,kB], intArr[1,kB], bd.wave[rL]] yArea = [ 0,0, a*intArr[1,kB]+b, bd.class[kC].spec[rL] ] ; yArea = [ 0,0, a*intArr[1,kB]+b, bd.spec[kC,rL] ] aR = poly_area( xArea,yArea) sR = intArr[1,kB]-bd.wave[rL] if debug eq 1 then cgPlot, /over, xArea, yArea, thick=2, color='orange' ;-possible place in the middle sM=0. aM=0. if lU lt rL then begin nMiss= rL-lU ; print, nMiss for kM=0, nMiss-1 do begin indX= lU+kM+[0,1,1,0] indY= lU+kM+[1,0] xArea = bd.wave[[indX]] yArea = [0,0,reform(bd.class[kC].spec[indY])] ; yArea = [0,0,reform(bd.spec[kC,indY])] aM += poly_area(xArea,yArea) sM += 0.050 if debug eq 1 then cgColorFill, xArea, yArea, color='gray' endfor endif area = aL + aR + aM length=sL + sR + sM refl[kC,kB] = area / length if debug eq 1 then cgPlot, /over, midW[kB], refl[kC,kB], psym='filled circle', symsize=3 endelse if ~finite(refl[kC,kB]) then begin print, '' print, bd.class[kC], track[kC,kB] print, intArr[*,kB] print, 'area : ', aL, aM, aR, area print, 'length: ', sL, sM, sR, length print, refl[kC,kB] stop endif ; kc=nbClass+1 endfor endfor refl[*,nbBand-1] = refl[*,nbBand-2] refl[*,0] = refl[*,1] ; for kC=0, nbClass-1 do begin ; ; cgPlot, bd.wave, bd.spec[kC,*], thick=4, title=bd.class[kC] ; ; cgPlot, vBand, [0,0]+0.15, thick=10, /over ; ; ; case1 = where( track[kC,*] eq 1 ) ; case2 = where( track[kC,*] eq 2 ) ; ; cgPlot, midW[case1], refl[kC,case1], color='green', thick=2, psym='filled circle', /over ; cgPlot, midW[case2], refl[kC,case2], color='blue', thick=2, psym='filled circle', /over ; ; ; endfor ;wdelete midW[0]=0.2 midW[nbBand-1]=250 openW, id, root+'lookup.dat', /get_lun line = '#waveIn waveMid waveOut Energy ' for kC=0, nbClass-1 do $ line += ' '+string(bd.class[kC],format='(A-4)') printf, id, line for kB=0, nbBand-1 do begin ; forprint, intArr[0,*], intArr[1,*], vArr line = string(intArr[0,kB],format='(F8.4)')+' '+$ string( midW[kB],format='(F8.4)')+' '+$ string(intArr[1,kB],format='(F8.4)')+' '+$ string(vArr[kB],format='(F7.4)')+' ' for kC=0, nbClass-1 do $ line += string( refl[kC,kB], format='(F6.4)')+' ' printf, id, line endfor close, id free_lun, id print, total(vArr) end