root = '/data/colors/sw-family/' dirDATA = root+'data/' xRang = [-.5,3] yRang = [-0.75, 0.25] nbClass=11 classDef = replicate({name:'',color:'',symbol:'',thick:1., size:1.}, nbClass) classDef.name = ['A','B','C','D','K','L','Q','R','S','V','X'] classDef.color= ['Brown', $ ; A 'Blue', $ ; B 'Black', $ ; C 'Gold', $ ; D 'Limegreen', $; K 'Skyblue',$ ; L 'Goldenrod',$ ; Q 'Purple', $ ; R 'Orange',$ ; S 'Red',$ ; V 'Gray'] ; X classDef.symbol=['Filled Up Triangle', $ ; A 'Filled Square', $ ; B 'Filled Circle', $ ; C 'Filled Down Triangle',$ ; D 'Filled Square', $ ; K 'Open Square', $ ; L 'Open Diamond', $ ; Q 'Open Down Triangle', $ ; R 'Star', $ ; S 'Open Right Triangle', $ ; V 'Filled Left triangle'] ; X classDef[8].size=1.5 ; C B X D L S K R Q V A slMin=[-0.50,-0.5, 0.2 , 0.7 , 0.7 , 0.6 , 0.6 , 0.5 , 0.5 , 0.6 , 2.15] slMax=[ 0.60, 0.0, 0.9 , 2.3 , 2.3 , 2.15, 1.2 , 2.3 , 0.95, 2.3 , 2.8 ] ziMin=[-0.20,-0.2,-0.03, 0.03 ,-0.1 ,-0.28, -0.11, -0.35,-0.28,-0.7 ,-0.3 ] ziMax=[ 0.15, 0.0, 0.15, 0.215, 0.03,-0.05, -0.04, -0.28,-0.2 ,-0.35,-0.1 ] ;shifting z-i up to account for difference between transform and no transform for loop=0,n_elements(ziMin)-1 do begin if loop eq 1 then goto, skipB ziMin(loop)=ziMin(loop)+0.035 ziMax(loop)=ziMax(loop)+0.035 skipB: endfor ord = [2,1,10,3,5,8,4,7,6,9,0] readcol, dirDATA+'all_families.csv', delimiter=',', /Silent, $ des,IAUNumber,IAUName,$ gMag,gUnc,rMag,rUnc,iMag,iUnc,zMag,zUnc,Phase,$ AlbedoV,SigAlbedoV,AlbedoNIR,SigAlbedoNIR,Diameter,SigDiameter,$ AbsoluteMagnitude,SemimajorAxis,InclinationEcliptic,OrbitalEccentricity,$ pSemimajorAxis,pInclination,pSinInclination,pOrbitalEccentricity, $ famId, famName, $ format='(A,L,A,'+$ 'F,F,F,F,F,F,F,F,F,'+$ 'F,F,F,F,F,F,'+$ 'F,D,F,F,'+$ 'D,F,F,F,L,A)' nbSSO = n_elements(IAUName) listFam = uniq(famId,sort(famId)) uniqFam=famId[listFam] nbFam=n_elements(listFam) ; values from Holmberg et al 2006 (MNRAS) ;-Solar Colours SunUmG = 1.40 & SunUmG_err = 0.08 SunGmG = 0.00 & SunGmG_err = 0.00 SunRmG = -0.45 & SunRmG_err = 0.02 SunImG = -0.55 & SunImG_err = 0.03 SunZmG = -0.61 & SunZmG_err = 0.04 wave=[0.4686, 0.6166, 0.7480, 0.8932 ] ;--III.1-- Convert UGRIZ Colors into Reflectance refG = 10.^( -0.4* ((gMag-gMag)-SunGmG) ) refR = 10.^( -0.4* ((rMag-gMag)-SunRmG) ) refI = 10.^( -0.4* ((iMag-gMag)-SunImG) ) refZ = 10.^( -0.4* ((zMag-gMag)-SunZmG) ) errG = abs(0.4*alog(10.)*refG * (gUnc + SunGmG_err) ) errR = abs(0.4*alog(10.)*refR * (rUnc + SunRmG_err) ) errI = abs(0.4*alog(10.)*refI * (iUnc + SunImG_err) ) errZ = abs(0.4*alog(10.)*refZ * (zUnc + SunZmG_err) ) gri = replicate({val:0.,unc:0.}, nbSSO) zi = replicate({val:0.,unc:0.}, nbSSO) class = replicate('U',nbSSO) ; window, 2 for kS=0, nbSSO-1 do begin if gMag[kS] lt 30 and $ rMag[kS] lt 30 and $ iMag[kS] lt 30 and $ zMag[kS] lt 30 then begin refl = [refG[kS], refR[kS], refI[kS], refZ[kS]] unc = [errG[kS], errR[kS], errI[kS], errZ[kS]] lin3 = linfit( wave[0:2],refl[0:2], measure_errors=unc[0:2], sigma=sig3 ) sig3[1]/=3. ; ; ww = [0.4, 0.8] ; lMid = lin3[0] + ww*lin3[1] ; lTop = lin3[0] + ww*(lin3[1]+sig3[1]) ; lBot = lin3[0] + ww*(lin3[1]-sig3[1]) ; ; xWW = [ww, reverse(ww)] ; yLL = [lTop, reverse(lBot)] ; ; cgPlot, wave, refl, yRange=[0.5,1.8], $ ; err_yLow=unc, err_yHigh=unc ; ; cgColorFill, xWW, yLL, color='light gray' ; cgPlot, /OverPlot, ww, lTop, color='Dark Gray' ; cgPlot, /OverPlot, ww, lBot, color='Dark Gray' ; cgPlot, /OverPlot, ww, lMid, color='CornFlower Blue' ; ; cgPlot, /OverPlot, wave, refl, err_yLow=unc, err_yHigh=unc gri[kS].val = lin3[1] gri[kS].unc = sig3[1] zi[kS].val = refZ[ks]-refI[kS] zi[kS].unc = errZ[ks]-errI[kS] for kC=0, nbClass-1 do begin if (gri[kS].val ge slMin[kC]) and $ (gri[kS].val lt slMax[kC]) and $ (zi[kS].val ge ziMin[kC]) and $ (zi[kS].val lt ziMax[kC]) then begin class[kS] = classDef[ord[kC]].name ; default is "U" if not in any class endif endfor endif endfor DIM = 1000. classArea = fltarr(DIM,DIM) xLen = xRang[1]-xRang[0] yLen = yRang[1]-yRang[0] for kC=0, nbClass-1 do begin xArr = [slMin[kC], slMax[kC], slMax[kC], slMin[kC], slMin[kC]] yArr = [ziMin[kC], ziMin[kC], ziMax[kC], ziMax[kC], ziMin[kC]] area = polyfillv( DIM*(xArr-xRang[0])/xLen, $ DIM*(yArr-yRang[0])/yLen, DIM, DIM) classArea[area] = kC+1 endfor cgPS_open, Filename=dirDATA+'fig-taxomembers.eps', $ /metric, /decomposed, /portrait, $ xsize=28, ysize=20, language_level=2, /quiet pos = [3500,2000,27500,19800] cgPlot, gri.val, zi.val, xRange=xRang, yRange=yRang, $ xTitle='gri slope (%/100 nm)', yTitle='z-i', $ pSym='Filled Circle', color='Light Gray', $ /Device, position=pos, /NoData ;--- three level of gray cgPlot, gri.val, zi.val, pSym='Filled Circle', /OverPlot, color='Light Gray' good = where( gri.unc le 3 and zi.unc le 0.03 ) cgPlot, gri[good].val, zi[good].val, /OverPlot, pSym='Filled Circle', color='Dark Gray' good = where( gri.unc le 1 and zi.unc le 0.01 ) cgPlot, gri[good].val, zi[good].val, /OverPlot, pSym='Filled Circle', color='Black' ;-- colors for classes totClassed=0 for kC=0, nbClass-1 do begin loc = where( strCmp(class,classDef[kC].name), tt ) print, classDef[kC].name, tt totClassed+=tt ; cgPlot, gri[loc].val, zi[loc].val, /OverPlot, $ ; pSym=classDef[kC].symbol, color=classDef[kC].color, symSize=classDef[kC].size endfor print, ' ',totClassed, nbSSO ; famCol=['Red','Green','Blue', $ ; 'Orange', 'Cornflower blue', 'YGB5', $ ; 'Chocolate', 'Yellow', 'Dark gray', $ ; 'PUR5', 'Navy', 'Beige', $ ; 'RED4', 'Magenta', 'RED7', $ ; 'Lawn Green', 'Thistle', 'Black', $ ; 'TG3'] ; for kF=0, nbFam-1 do begin ; loc = where( famId eq uniqFam[kF] ) ; cgPlot, gri[loc].val, zi[loc].val, /OverPlot, $ ; pSym='Filled Circle', color=famCol[kF] ; ; ; xKey = 2.5 ; yKey= -0.7 ; yShi = 0.05 ; ; cgPlot, /OverPlot, xKey, yKey+kF*yShi, pSym='Filled Circle', color=famCol[kF] ; cgText, xKey+0.1, yKey+kF*yShi-0.005, famName[listFam[kF]] ; ; ; endfor ;--- draw zones for kC=0, nbClass-1 do begin classZone = fltarr(DIM,DIM) area = where( classArea eq kC+1 ) classZone[area] = kC+1 contour, classZone, $ min_value=0, max_value=nbClass+1, $ levels=[kC+1], $ /Closed, /Path_Data_Coords, $ /Device, /Noerase, $ xSt=1, ySt=1, $ position = pos, $ path_info=INFO, path_xy=XY S = [indgen(INFO[0].N), 0] coordXY=transpose(xy[*,INFO[0].OFFSET + S ]) cgPlot, coordXY[*,0], coordXY[*,1], color=classDef[ord[kC]].color, thick=2, /OverPlot cgText, max(coordXY[*,0])-30, min(coordXY[*,1])+10, $ classDef[ord[kC]].name, $ charsize=1.5, charth=2, color=classDef[ord[kC]].color endfor cgPS_close, /PNG ;do some selection here: value and uncertainties below threshold X select = where( gri.unc le 1 and zi.unc le 0.01 and $ gri.val le 3.5 and gri.val ge -1 and $ zi.val ge -1 and zi.val le 0.5 and $ gMag lt 30 and $ rMag lt 30 and $ iMag lt 30 and $ zMag lt 30, nbSel ) print, nbSel forprint, des,class, gri.val, gri.unc, zi.val, zi.unc, $ textout=dirDATA+'all_families-taxoslope.csv', /Silent, $ format='(A,",",A,",",F,",",F,",",F,",",F)', $ subset=select, $ comment='Designation,Class, slope, dSlope, ZI, dZI' end