root = '/home/bcarry/work/data/q-s-types/' file = 'svo-uniq.dat' cwd=root+'round1/' ; file = 'mba-uniq.dat' ; cwd=root+'round2/' doF1 = 01 ;-- Perihelion vs Slope doF2 = 0 ;-- Perihelion vs Fraction doF3 = 0 ;-- Average vs Slope doF4 = 0 ;-- Average vs Fraction doF5 = 01 ;-- H distribution & H vs Fraction doF6 = 0 ;-- Flux vs Slope doF7 = 01 ;-- Slope vs H svo = sdssReadTaxonomy( root+file ) peri = svo.a*( 1-svo.e ) avg = svo.a*( 1 + svo.e*svo.e/2. ) flux = 1/(svo.a*svo.a*sqrt( 1-svo.e*svo.e) ) nbSVO = n_elements( svo ) q=where( strCmp( svo.class, 'Q', 1) and svo.taxo.sl gt 0 and svo.taxo.sl le 2.5, nbQ ) s=where( strCmp( svo.class, 'S', 1) and svo.taxo.sl gt 0 and svo.taxo.sl le 2.5, nbS ) b=[s,q] print, nbQ, nbS, nbSVO aArr = [0.3,2] ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------First figure: perihelion vs slope if doF1 eq 1 then begin print, '' print, ' --1-- Perihelion vs Slope' cgPS_open, Filename=cwd+'qs-peri_vs_slope.eps', /metric, /decomposed, /encapsulated, $ xsize=30, ysize=15, language_level=2, /quiet cgPlot, peri[s], svo[s].taxo.sl, /NoData, $ xStyle=1, xRange=[0.2,1.8], $ yStyle=1, yRange=[0,2.5], $ xTitle='Perihelion (au)', $ yTitle='Spectral slope' cgPlot, /OverPlot, 0.39*[1,1], [0,100], color='Gray', lineStyle=2 cgPlot, /OverPlot, 0.72*[1,1], [0,100], color='Gray', lineStyle=2 cgPlot, /OverPlot, 1.00*[1,1], [0,100], color='Gray', lineStyle=2 cgPlot, /OverPlot, 1.53*[1,1], [0,100], color='Gray', lineStyle=2 cgText, 0.39+0.01, 2.2, color='Gray', 'M' cgText, 0.72+0.01, 2.2, color='Gray', 'V' cgText, 1.00+0.01, 2.2, color='Gray', 'E' cgText, 1.53+0.01, 2.2, color='Gray', 'M' cgPlot, /OverPlot, peri[s], svo[s].taxo.sl, psym='Filled circle' aS = regress( peri[s], svo[s].taxo.sl, const=bS, yFit=sFit ) cgPlot, /OverPlot, peri[s], sFit print, '--- S-types ---' print, 'S0: ', bS print, 'Ss: ', aS cgPlot, /OverPlot, peri[q], svo[q].taxo.sl, psym='Filled square', color='Red', symSize=1.5 aQ = regress( reform(peri[q]), reform(svo[q].taxo.sl), const=bQ, yFit=qFit ) cgPlot, /OverPlot, peri[q], qFit, color='Red' print, '--- Q-types ---' print, 'S0: ', bQ print, 'Sq: ', aQ ; cgPlot, /OverPlot, peri[b], svo[b].taxo.sl, psym=' circle' ; aB = regress( peri[b], svo[b].taxo.sl, const=bB, yFit=sFit ) ; cgPlot, /OverPlot, peri[b], sFit, color='Blue' ; ; print, '--- Both types ---' ; print, 'S0: ', bB ; print, 'Sb: ', aB low = where( peri[b] le 1., comp=high ) aBL = regress( peri[b[low]], svo[b[low]].taxo.sl, const=bBL, yFit=sFit, correl=cBL ) cgPlot, /OverPlot, peri[b[low]], sFit, color='Blue' print, '--- Both Low ---' print, 'S0: ', bBL print, 'Sb: ', aBL print, 'c : ', cBL aBH = regress( peri[b[high]], svo[b[high]].taxo.sl, const=bBH, yFit=sFit, correl=cBH ) cgPlot, /OverPlot, peri[b[high]], sFit, color='Blue' print, '--- Both High ---' print, 'S0: ', bBH print, 'Sb: ', aBH print, 'c : ', cBH cgPS_close, /png, Delete_PS=0 endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------Second figure: Q/S ratio if doF2 eq 1 then begin print, '' print, ' --2-- Perihelion vs Fraction' cgHistoPlot, peri[q], binSize=0.1, loc=qArr, histData=histoQ, minInput=0.3, maxInput=1.9 cgPlot, /OverPlot, qArr, histoQ, psym=4 ; stop cgHistoPlot, peri[s], binSize=0.1, loc=sArr, histData=histoS, minInput=0.3, maxInput=1.9 cgPlot, /OverPlot, sArr, histoS, psym=4 cgPlot, sArr, qArr ; stop wdelete cgPS_open, Filename=cwd+'qs-peri_vs_fraction.eps', /metric, /decomposed, /encapsulated, $ xsize=30, ysize=15, language_level=2, /quiet cgPlot, qArr, 100.*(histoQ)/histoS, /NoData, $ xStyle=1, xRange=[0.2,1.8], $ yStyle=1, yRange=[0,100], $ xTitle='Perihelion (au)', $ yTitle='Q/S fraction (%)' cgPlot, /OverPlot, 0.39*[1,1], [0,100], color='Gray', lineStyle=2 cgPlot, /OverPlot, 0.72*[1,1], [0,100], color='Gray', lineStyle=2 cgPlot, /OverPlot, 1.00*[1,1], [0,100], color='Gray', lineStyle=2 cgPlot, /OverPlot, 1.53*[1,1], [0,100], color='Gray', lineStyle=2 cgText, 0.39+0.01, 90, color='Gray', 'M' cgText, 0.72+0.01, 90, color='Gray', 'V' cgText, 1.00+0.01, 90, color='Gray', 'E' cgText, 1.53+0.01, 90, color='Gray', 'M' ; cgPlot, /OverPlot, qArr, 800.*histoS/max(histoS), psym=10, color='Light gray' cgPlot, /OverPlot, qArr, 100.*histoQ/histoS, psym=10 cgPS_close, /png, Delete_PS=0 endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------Third figure: average vs slope if doF3 eq 1 then begin print, '' print, ' --3-- Average vs Slope' cgPS_open, Filename=cwd+'qs-avg_vs_slope.eps', /metric, /decomposed, /encapsulated, $ xsize=30, ysize=15, language_level=2, /quiet cgPlot, avg[s], svo[s].taxo.sl, /NoData, $ xStyle=1, xRange=[0.9,3.5], $ yStyle=1, yRange=[0,2.5], $ xTitle='Average distance (au)', $ yTitle='Spectral slope' cgPlot, /OverPlot, 0.39*[1,1], [0,100], color='Gray', lineStyle=2 cgPlot, /OverPlot, 0.72*[1,1], [0,100], color='Gray', lineStyle=2 cgPlot, /OverPlot, 1.00*[1,1], [0,100], color='Gray', lineStyle=2 cgPlot, /OverPlot, 1.53*[1,1], [0,100], color='Gray', lineStyle=2 cgText, 0.39+0.01, 2.2, color='Gray', 'M' cgText, 0.72+0.01, 2.2, color='Gray', 'V' cgText, 1.00+0.01, 2.2, color='Gray', 'E' cgText, 1.53+0.01, 2.2, color='Gray', 'M' cgPlot, /OverPlot, avg[s], svo[s].taxo.sl, psym='Filled circle' aS = regress( avg[s], svo[s].taxo.sl, const=bS, yFit=sFit ) cgPlot, /OverPlot, avg[s], sFit print, '--- S-types ---' print, 'S0: ', bS print, 'Ss: ', aS cgPlot, /OverPlot, avg[q], svo[q].taxo.sl, psym='Filled square', color='Red', symSize=1.5 aQ = regress( reform(avg[q]), reform(svo[q].taxo.sl), const=bQ, yFit=qFit ) cgPlot, /OverPlot, avg[q], qFit, color='Red' print, '--- Q-types ---' print, 'S0: ', bQ print, 'Sq: ', aQ cgPS_close, /png, Delete_PS=0 endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; if doF4 eq 1 then begin print, '' print, ' --4-- Average vs Fraction' cgHistoPlot, avg[q], binSize=0.1, loc=qArr, histData=histoQ, minInput=0.9, maxInput=3.9 cgPlot, /OverPlot, qArr, histoQ, psym=4 ; stop cgHistoPlot, avg[s], binSize=0.1, loc=sArr, histData=histoS, minInput=0.9, maxInput=3.9 cgPlot, /OverPlot, sArr, histoS, psym=4 cgPlot, sArr, qArr ; cgPlot, sArr, qArr ; stop wdelete cgPS_open, Filename=cwd+'qs-avg_vs_fraction.eps', /metric, /decomposed, /encapsulated, $ xsize=30, ysize=15, language_level=2, /quiet cgPlot, qArr, 100.*(histoQ)/histoS, /NoData, $ xStyle=1, xRange=[0.9,3.5], $ yStyle=1, yRange=[0,210], $ xTitle='Average distance (au)', $ yTitle='Q/S fraction (%)' cgPlot, /OverPlot, 0.39*[1,1], [0,1000], color='Gray', lineStyle=2 cgPlot, /OverPlot, 0.72*[1,1], [0,1000], color='Gray', lineStyle=2 cgPlot, /OverPlot, 1.00*[1,1], [0,1000], color='Gray', lineStyle=2 cgPlot, /OverPlot, 1.53*[1,1], [0,1000], color='Gray', lineStyle=2 ; cgText, 0.39+0.01, 185, color='Gray', 'M' ; cgText, 0.72+0.01, 185, color='Gray', 'V' cgText, 1.00+0.01, 185, color='Gray', 'E' cgText, 1.53+0.01, 185, color='Gray', 'M' cgPlot, /OverPlot, qArr, 100.*histoQ/histoS, psym=10 cgPS_close, /png, Delete_PS=0 endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; if doF5 eq 1 then begin print, '' print, ' --4-- Histo H and H vs Fraction' binS = 1 cgHistoPlot, svo[q].h, binSize=binS, loc=qArr, histData=histoQ, minInput=10, maxInput=25 cgPlot, /OverPlot, qArr, histoQ, psym=4 ; stop cgHistoPlot, svo[s].h, binSize=binS, loc=sArr, histData=histoS, minInput=10, maxInput=25 cgPlot, /OverPlot, sArr, histoS, psym=4 ; stop wdelete cgPS_open, Filename=cwd+'qs-H.eps', /metric, /decomposed, /encapsulated, $ xsize=30, ysize=15, language_level=2, /quiet cgPlot, sArr, 100.*histoS/total(histoS), psym=10, $ xTitle='Absolute magnitude H', $ yTitle='Fraction of population (%)' cgPlot, /OverPlot, qArr, 100.*histoQ/total(histoQ), psym=10, color='Red' cgPS_close, /png, Delete_PS=0 cgPS_open, Filename=cwd+'qs-H_vs_fraction.eps', /metric, /decomposed, /encapsulated, $ xsize=30, ysize=15, language_level=2, /quiet cgPlot, qArr, 100.*(histoQ)/histoS, /NoData, $ xStyle=1, xRange=[12,23], $ yStyle=1, yRange=[0,110], $ xTitle='Absolute magnitude H', $ yTitle='Q/S fraction (%)' cgPlot, /OverPlot, qArr, 100.*histoQ/histoS, psym=10 cgPS_close, /png, Delete_PS=0 endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; if doF6 eq 01 then begin print, '' print, ' --6-- Flux vs Slope' cgPS_open, Filename=cwd+'qs-flux_vs_slope.eps', /metric, /decomposed, /encapsulated, $ xsize=30, ysize=15, language_level=2, /quiet cgPlot, flux[s], svo[s].taxo.sl, /NoData, $ xStyle=1, xRange=[0.1,1.8], $ yStyle=1, yRange=[0,2.5], $ xTitle='Flux over an orbit', $ yTitle='Spectral slope' cgPlot, /OverPlot, flux[s], svo[s].taxo.sl, psym='Filled circle' aS = regress( flux[s], svo[s].taxo.sl, const=bS, yFit=sFit ) cgPlot, /OverPlot, flux[s], sFit print, '--- S-types ---' print, 'S0: ', bS print, 'Ss: ', aS cgPlot, /OverPlot, flux[q], svo[q].taxo.sl, psym='Filled square', color='Red', symSize=1.5 aQ = regress( reform(flux[q]), reform(svo[q].taxo.sl), const=bQ, yFit=qFit ) cgPlot, /OverPlot, flux[q], qFit, color='Red' print, '--- Q-types ---' print, 'S0: ', bQ print, 'Sq: ', aQ ; cgPlot, /OverPlot, flux[b], svo[b].taxo.sl, psym=' circle' aB = regress( flux[b], svo[b].taxo.sl, const=bB, yFit=sFit ) cgPlot, /OverPlot, flux[b], sFit, color='Blue' print, '--- Both types ---' print, 'S0: ', bB print, 'Sb: ', aB cgPS_close, /png, Delete_PS=0 endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; if doF7 eq 01 then begin print, '' print, ' --7-- H vs Slope' cgPS_open, Filename=cwd+'qs-H_vs_slope.eps', /metric, /decomposed, /encapsulated, $ xsize=30, ysize=15, language_level=2, /quiet cgPlot, svo[s].h, svo[s].taxo.sl, /NoData, $ xStyle=1, xRange=[10,25], $ yStyle=1, yRange=[0,2.5], $ xTitle='Absolute magnitude (H)', $ yTitle='Spectral slope' cgPlot, /OverPlot, svo[s].h, svo[s].taxo.sl, psym='Filled circle' aS = regress( svo[s].h, svo[s].taxo.sl, const=bS, yFit=sFit ) cgPlot, /OverPlot, svo[s].h, sFit, color='gray' print, '--- S-types ---' print, 'S0: ', bS print, 'Ss: ', aS cgPlot, /OverPlot, svo[q].h, svo[q].taxo.sl, psym='Filled square', color='Red', symSize=1.5 aQ = regress( reform(svo[q].h), reform(svo[q].taxo.sl), const=bQ, yFit=qFit ) cgPlot, /OverPlot, svo[q].h, qFit, color='Red' print, '--- Q-types ---' print, 'S0: ', bQ print, 'Sq: ', aQ ; cgPlot, /OverPlot, flux[b], svo[b].taxo.sl, psym=' circle' ; aB = regress( svo[b].h, svo[b].taxo.sl, const=bB, yFit=sFit, correl=cB ) ; cgPlot, /OverPlot, svo[b].h, sFit, color='Blue' ; ; print, '--- Both types ---' ; print, 'S0: ', bB ; print, 'Sb: ', aB ; print, 'C: ', cB cgPS_close, /png, Delete_PS=0 endif end