;--- Euclid G2V colors: ; V-Vis = 0.528522 ; V-Y = 0.954786 ; V-J = 1.28468 ; V-H = 1.48595 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Directories spawn, 'hostname', host case strTrim(host) of 'hyperion': dirEuclid = '/observ/euclid/' 'endymion': dirEuclid = '/home/bcarry/work/data/euclid/' else: stop endcase filterDir = dirEuclid+'filters/' dirTaxo = dirEuclid+'taxo/' ;-What to do showGraphBDM = 0 ;-X createGraphBDM = 01 ;-PNG/EPS showGraphMOVIS = 0 ;-X createGraphMOVIS = 01 ;-PNG/EPS checkTaxo = 01 ;-Classify BDM colors instead of MOVIS (-> auto-check) onlyValidFig = 0 ;-Re run barplot code only - only for checktaxo=1 ;-Plot Parameters ccRang = [ [ 0.10, 0.90], $ ;- Vis-Y [ 0.45, 1.40], $ ;- Vis-J [ 0.35, 1.40], $ ;- Vis-H ; [ 0.45, 1.60], $ ;- Vis-H [ 0.00, 0.90], $ ;- Y-J [-0.10, 1.10], $ ;- Y-H [-0.15, 0.30] ] ;- J-H ccTickInt= [0.2, 0.2, 0.2, 0.2, 0.2, 0.1 ] ccTickMin= [ 4 , 4 , 4 , 4 , 4 , 4 ] numSel=3 ; if onlyValidFig eq 1 then goto, j2xm ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- I --- Define classes ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--I.1--Class Labels -------------------------------------------------------------------- cLabel=['A','B','C','D','K','L','Q','S','T','X','V'] nbClass = n_elements( cLabel ) print, 'List of classes considered here ('+strtrim(string(nbClass),2)+') ' print, cLabel ;--I.2--Class Look'n'Feel --------------------------------------------------------------- cColor=['Sienna' , $ ;-- A 'Gray' , $ ;-- B 'Charcoal' , $ ;-- C 'Dark Gray' , $ ;-- D 'Green' , $ ;-- K 'Navy', $ ;-- L 'Orange' , $ ;-- Q 'Cornflower blue', $ ;-- S 'PUR8' , $ ;-- T 'Magenta' , $ ;-- X 'Red' ] ;-- V cSymb =['Filled circle' , $ ;-- A 'Filled circle', $ ;-- B 'Open circle', $ ;-- C 'Filled circle', $ ;-- D 'Open circle', $ ;-- K 'Open square', $ ;-- L 'Filled star', $ ;-- Q 'Open square', $ ;-- S 'Star' , $ ;-- T 'Open diamond', $ ;-- X 'Filled Circle' ] ;-- V SSymb =[ 0.8 ,$ ;-- A 0.8 ,$ ;-- B 0.8 ,$ ;-- C 0.8 ,$ ;-- D 0.8 ,$ ;-- K 0.8 ,$ ;-- L 1.5 ,$ ;-- Q 0.8 ,$ ;-- S 1.5 ,$ ;-- T 0.8 ,$ ;-- X 0.8 ] ;-- V ;--I.X--Useful stuff -------------------------------------------------------------------- theta=2*!PI*findgen(361)/360 ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- II --- Read Datasets ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--II.1-- Solar Color in Euclid [from ESA sky stupid filter by hand by me] ;--- Euclid G2V colors: ; V-Vis = 0.528522 ; V-Y = 0.954786 ; V-J = 1.28468 ; V-H = 1.48595 sunVisY = 0.954786 - 0.528522 sunVisJ = 1.28468 - 0.528522 sunVisH = 1.48595 - 0.528522 sunY__J = 1.28468 - 0.954786 sunY__H = 1.48595 - 0.954786 sunJ__H = 1.48595 - 1.28468 print, '' print, ' II. Solar colors in Euclid' print, ' sunVisY', sunVisY print, ' sunVisJ', sunVisJ print, ' sunVisH', sunVisH print, ' sunY__J', sunY__J print, ' sunY__H', sunY__H print, ' sunJ__H', sunJ__H ;--II.2--Read Bus-DeMeo Summary Color --------------------------------------------------- ;--II.2.1--Read Bus-DeMeo reflectance file readcol, dirTaxo+'euclid-demeoColors.csv', $ ssoNum, ssoName, ssoClass, $ colVis, colY, colJ, colH, $ uncVis, uncY, uncJ, uncH, $ classId, compLabel, complexId, $ format='(F6,A-15,A-3,F7.4,F7.4,F7.4,F7.4,F7.4,F7.4,F7.4,F7.4,I2,A1,I1)', $ skipline=1, delim=',', /Silen nbSSO= n_elements( ssoNum ) ;--II.2.2--Convert to magnitude VmY = 2.5*alog10( colY ) + sunVisY VmJ = 2.5*alog10( colJ ) + sunVisJ VmH = 2.5*alog10( colH ) + sunVisH YmJ = 2.5*alog10( colJ/colY ) + sunY__J YmH = 2.5*alog10( colH/colY ) + sunY__H JmH = 2.5*alog10( colH/colJ ) + sunJ__H ;--II.2.3--Uncertainties dVmY = abs(2.5*alog10( colY+uncY ) - 2.5*alog10( colY-uncY )) dVmJ = abs(2.5*alog10( colJ+uncJ ) - 2.5*alog10( colJ-uncJ )) dVmH = abs(2.5*alog10( colH+uncH ) - 2.5*alog10( colH-uncH )) dYmJ = abs(2.5*alog10( (colJ+uncJ)/(colY-uncY) ) - 2.5*alog10( (colJ-uncJ)/(colY+uncY) )) dYmH = abs(2.5*alog10( (colH+uncH)/(colY-uncY) ) - 2.5*alog10( (colH-uncH)/(colY+uncY) )) dJmH = abs(2.5*alog10( (colH+uncH)/(colJ-uncJ) ) - 2.5*alog10( (colH-uncH)/(colJ+uncJ) )) ;--II.2.4-- Write down the magnitudes in CSV forprint, ssoNum, ssoName, ssoClass, $ VmY, VmJ, VmH, YmJ, YmH, JmH, $ dVmY, dVmJ, dVmH, dYmJ, dYmH, dJmH, $ classId, compLabel, complexId, $ format='(I-6,",",A-15,",",A-3,12(",",F7.4),",",I2,",",A1,",",I2)', $ textout=dirTaxo+'euclid-ssomag.csv', $ comment='Num, Name, Class, '+$ '(VIS-Y), (VIS-J), (VIS-H), (Y-J), (Y-H), (J-H), '+$ 'd(VIS-Y), d(VIS-J), d(VIS-H), d(Y-J), d(Y-H), d(J-H), '+$ 'classId,complex,complexId', $ /silent ;--II.2.5-- Analyze # of classes bdmCol=[[VmY],[VmJ],[VmH],[YmJ],[YmH],[JmH]] class = compLabel uClass = uniq(class, sort(class)) print, 'List of classes in Bus-DeMeo color file ('+strtrim(string(n_elements(uClass)),2)+') ' print, class[uClass] ;--II.3--Euclid Filters 7 Combinations --------------------------------------------- filters=['VIS', 'Y', 'J', 'H'] nbFilt=n_elements(filters) nbComb = factorial(nbFilt)/ ( 2*factorial( nbFilt-2 )) label=['VIS-Y','VIS-J','VIS-H','Y-J','Y-H','J-H'] ;---------- TEST ;--II.2.5-- Analyze # of classes bdmCol=[[VmY],[VmJ],[YmH]] class = compLabel uClass = uniq(class, sort(class)) print, 'List of classes in Bus-DeMeo color file ('+strtrim(string(n_elements(uClass)),2)+') ' print, class[uClass] ;--II.3--Euclid Filters 7 Combinations --------------------------------------------- filters=['VIS', 'Y', 'J', 'H'] nbFilt=n_elements(filters) nbComb = 3 label=['VIS-Y','VIS-J','Y-H'] ;--II.4--Autocheck ---------------------------------------------------------------------- ssoCol = bdmCol nbSSO = n_elements( class ) print, 'Number of asteroids in Bus-DeMeo:', nbSSO ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- III --- Compute Mean + StdDev for Each Class ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--III.1--Output Structures ------------------------------------------------------------- covArr = replicate({correl:0., covar:0., angle:0., cos:0., sin:0., $ meanX:0., meanY:0., stdX:0., stdY:0., stdA:0., stdB:0.} ,nbClass,nbComb,nbComb) for kClass=0, nbClass-1 do begin ;--III.2--Select Only Asteroids from Current Class ----------------------------------- placeCl = where( strCmp(class, cLabel[kClass]), nbCl) print, cLabel[kClass]+'-types ('+strtrim(string(nbCl),2)+')' ;--III.3--Compute Correlations and Covariances --------------------------------------- locBDM = transpose(bdmCol[ placeCl,* ]) corr = correlate(locBDM) cov = correlate(locBDM, /cov ) ;--III.3--Determine Could Orientation and Deviation ---------------------------------- for kComb1=0, nbComb-1 do begin ;--III.3.1--Sample Variance & Normalization of Covariance varComb1 = variance( bdmCol[placeCl,kComb1] ) cov[kComb1,*] /= varComb1 ;--III.3.2--Loop over Second Dimension for kComb2=kComb1+1, nbComb-1 do begin ;--III.3.3--Keep Correlation and Covariance covArr[kClass,kComb1,kComb2].correl = corr[kComb1,kComb2] covArr[kClass,kComb1,kComb2].covar = cov[kComb1,kComb2] ;--III.3.4--Compute Cloud Center covArr[kClass,kComb1,kComb2].meanX = mean( bdmCol[ placeCl,kComb1] ) covArr[kClass,kComb1,kComb2].meanY = mean( bdmCol[ placeCl,kComb2] ) covArr[kClass,kComb1,kComb2].stdX = stddev( bdmCol[ placeCl,kComb1] ) covArr[kClass,kComb1,kComb2].stdY = stddev( bdmCol[ placeCl,kComb2] ) ;--III.3.5--Compute Cloud Orientation covArr[kClass,kComb1,kComb2].angle = atan( covArr[kClass,kComb1,kComb2].covar ) covArr[kClass,kComb1,kComb2].cos = cos( covArr[kClass,kComb1,kComb2].angle ) covArr[kClass,kComb1,kComb2].sin = sin( covArr[kClass,kComb1,kComb2].angle ) ;--III.3.6--Compute Cloud Dispersion newX = bdmCol[placeCl,kComb1] - covArr[kClass,kComb1,kComb2].meanX newY = bdmCol[placeCl,kComb2] - covArr[kClass,kComb1,kComb2].meanY rotX = newX*covArr[kClass,kComb1,kComb2].cos + newY*covArr[kClass,kComb1,kComb2].sin rotY =-newX*covArr[kClass,kComb1,kComb2].sin + newY*covArr[kClass,kComb1,kComb2].cos covArr[kClass,kComb1,kComb2].stdA = stddev( rotX ) covArr[kClass,kComb1,kComb2].stdB = stddev( rotY ) endfor endfor endfor ;--III.3--Validation Figures on Screen --------------------------------------------------- if showGraphBDM eq 1 then begin ;--III.4.1-- Prepare plot cgPS_open, Filename=dirTaxo+'DeMeo-Taxonomy.eps', /metric, /decomposed, /encapsulated, $ xSize=30, ySize=30, language_level=2, /quiet xLen = 0.44 ; 0.176 yLen = xLen xSep = 0.01 ySep = 0.01 xBot = 0.085 yBot = 0.085 ;--III.4.2-- Loop over Combinations for kComb1=0, nbComb-1 do begin for kComb2=nbComb-1, kComb1+1, -1 do begin ;--III.4.3--Prepare Plot OutLay ;--III.4.3/A--XY Range lonX = [ covArr[*,kComb1,kComb2].meanX+covArr[*,kComb1,kComb2].stdA, $ covArr[*,kComb1,kComb2].meanX-covArr[*,kComb1,kComb2].stdA, bdmCol[*,kComb1]] lonY = [ covArr[*,kComb1,kComb2].meanY+covArr[*,kComb1,kComb2].stdB, $ covArr[*,kComb1,kComb2].meanY-covArr[*,kComb1,kComb2].stdB, bdmCol[*,kComb2]] xRang = ccRang[*,kComb1] yRang = ccRang[*,kComb2] ;--III.4.3/B--Graph position pos= [xBot + kComb1*xLen, $ yBot + (kComb2-1)*xLen, $ xBot + (kComb1+1)*xLen, $ yBot + (kComb2 )*xLen ] ;--III.4.3/C--XY Axes On the First column if kComb1 eq 0 then begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTitle=label[kComb2], $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTitle=label[kComb2], $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] endelse ;--III.4.3/D--XY Axes On the other columns endif else begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20), $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20), $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] endelse endelse ;--III.4.4--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(class,cLabel[kClass]) ) cgPlot, /OverPlot, bdmCol[cur,kComb1], bdmCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--III.4.5--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin ;--III.4.5/A--Long-Axis Line xA = [-1,1]*covArr[kClass,kComb1,kComb2].stdA yA = [0,0] rotXA = xA*covArr[kClass,kComb1,kComb2].cos - yA*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYA = xA*covArr[kClass,kComb1,kComb2].sin + yA*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXA, rotYA, color=cColor[kClass], /OverPlot ;--III.4.5/B--Short-Axis Line xB = [0,0] yB = [-1,1]*covArr[kClass,kComb1,kComb2].stdB rotXB = xB*covArr[kClass,kComb1,kComb2].cos - yB*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYB = xB*covArr[kClass,kComb1,kComb2].sin + yB*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXB, rotYB, color=cColor[kClass], /OverPlot ;--III.4.5/C--Ellipse x = covArr[kClass,kComb1,kComb2].stdA * cos(theta) y = covArr[kClass,kComb1,kComb2].stdB * sin(theta) rotX = x*covArr[kClass,kComb1,kComb2].cos - y*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotY = x*covArr[kClass,kComb1,kComb2].sin + y*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotX, rotY, color=cColor[kClass], /OverPlot cgText, covArr[kClass,kComb1,kComb2].meanX+0.015, $ covArr[kClass,kComb1,kComb2].meanY+0.015, cLabel[kClass], color=cColor[kClass] endfor endfor endfor ;--III.4.5-- Display Key pos= [xBot + 1*xLen, $ 0 * yLen, $ xBot + 2*xLen, $ 1 * yLen ] cgPlot, [0,1], [0,1], /NoData, /NoErase, /Normal, position=pos, xStyle=5, yStyle=5 xKey = 0.2 xShi = 0.03 yKey = 0.90 yShi =-0.075 yTxt = -0.015 for kClass=0, nbClass-1 do begin cur=where( strcmp(class,cLabel[kClass]), nbCur ) cgPlot, xKey, yKey+kClass*yShi, psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass], /OverPlot cgText, xKey+xShi, yKey+kClass*yShi+yTxt, cLabel[kClass]+'-types: ', color=cColor[kClass] cgText, xKey+xShi*13, yKey+kClass*yShi+yTxt, strtrim(string(nbCur),2), align=1, color=cColor[kClass] endfor cgPS_close, /png, Delete_PS=0 ;--III.3.1--Loop over Filter Combinations for kComb1=0, nbComb-1 do begin for kComb2=kComb1+1, nbComb-1 do begin ;--III.3.2--Prepare Plot OutLay lonX = [ covArr[*,kComb1,kComb2].meanX+covArr[*,kComb1,kComb2].stdA, $ covArr[*,kComb1,kComb2].meanX-covArr[*,kComb1,kComb2].stdA, bdmCol[*,kComb1]] lonY = [ covArr[*,kComb1,kComb2].meanY+covArr[*,kComb1,kComb2].stdB, $ covArr[*,kComb1,kComb2].meanY-covArr[*,kComb1,kComb2].stdB, bdmCol[*,kComb2]] cgPlot, lonX, lonY, /NoData, xTitle=label[kComb1], yTitle=label[kComb2] ;--III.3.3--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(class,cLabel[kClass]) ) cgPlot, /OverPlot, bdmCol[cur,kComb1], bdmCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--III.3.4--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin cgPlot, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*cos(theta), $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*sin(theta), color=cColor[kClass], /over cgErrPlot, covArr[kClass,kComb1,kComb2].meanX, color=cColor[kClass], $ covArr[kClass,kComb1,kComb2].meanY-covArr[kClass,kComb1,kComb2].stdB, $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB cgErrPlot, covArr[kClass,kComb1,kComb2].meanY, color=cColor[kClass], $ covArr[kClass,kComb1,kComb2].meanX-covArr[kClass,kComb1,kComb2].stdA, $ covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA, /horizontal cgText, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*0.3, $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*0.3, cLabel[kClass], color=cColor[kClass] endfor ; stop endfor endfor endif ;--III.4--Validation Figures in a File -------------------------------------------------- if createGraphBDM eq 1 then begin ;--III.4.1-- Prepare plot cgPS_open, Filename=dirTaxo+'DeMeo-Taxonomy.eps', /metric, /decomposed, /encapsulated, $ xSize=40, ySize=40, language_level=2, /quiet xLen = 0.176 yLen = xLen xSep = 0.01 ySep = 0.01 xBot = 0.085 yBot = 0.05 ;--III.4.2-- Loop over Combinations for kComb1=0, nbComb-1 do begin for kComb2=nbComb-1, kComb1+1, -1 do begin ;--III.4.3--Prepare Plot OutLay ;--III.4.3/A--XY Range lonX = [ covArr[*,kComb1,kComb2].meanX+covArr[*,kComb1,kComb2].stdA, $ covArr[*,kComb1,kComb2].meanX-covArr[*,kComb1,kComb2].stdA, bdmCol[*,kComb1]] lonY = [ covArr[*,kComb1,kComb2].meanY+covArr[*,kComb1,kComb2].stdB, $ covArr[*,kComb1,kComb2].meanY-covArr[*,kComb1,kComb2].stdB, bdmCol[*,kComb2]] xRang = ccRang[*,kComb1] yRang = ccRang[*,kComb2] ;--III.4.3/B--Graph position pos= [xBot + kComb1*xLen, $ yBot + (kComb2-1)*xLen, $ xBot + (kComb1+1)*xLen, $ yBot + (kComb2 )*xLen ] ;--III.4.3/C--XY Axes On the First column if kComb1 eq 0 then begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTitle=label[kComb2], $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTitle=label[kComb2], $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] endelse ;--III.4.3/D--XY Axes On the other columns endif else begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20), $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20), $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] endelse endelse ;--III.4.4--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(class,cLabel[kClass]) ) cgPlot, /OverPlot, bdmCol[cur,kComb1], bdmCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--III.4.5--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin ;--III.4.5/A--Long-Axis Line xA = [-1,1]*covArr[kClass,kComb1,kComb2].stdA yA = [0,0] rotXA = xA*covArr[kClass,kComb1,kComb2].cos - yA*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYA = xA*covArr[kClass,kComb1,kComb2].sin + yA*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXA, rotYA, color=cColor[kClass], /OverPlot ;--III.4.5/B--Short-Axis Line xB = [0,0] yB = [-1,1]*covArr[kClass,kComb1,kComb2].stdB rotXB = xB*covArr[kClass,kComb1,kComb2].cos - yB*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYB = xB*covArr[kClass,kComb1,kComb2].sin + yB*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXB, rotYB, color=cColor[kClass], /OverPlot ;--III.4.5/C--Ellipse x = covArr[kClass,kComb1,kComb2].stdA * cos(theta) y = covArr[kClass,kComb1,kComb2].stdB * sin(theta) rotX = x*covArr[kClass,kComb1,kComb2].cos - y*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotY = x*covArr[kClass,kComb1,kComb2].sin + y*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotX, rotY, color=cColor[kClass], /OverPlot endfor endfor endfor ;--III.4.5-- Display Key pos= [xBot + 3*xLen, $ yBot + (0)*yLen, $ xBot + (5)*xLen, $ yBot + (2 )*yLen ] cgPlot, [0,1], [0,1], /NoData, /NoErase, /Normal, position=pos, xStyle=5, yStyle=5 xKey = 0.2 xShi = 0.03 yKey = 0.90 yShi =-0.075 yTxt = -0.015 for kClass=0, nbClass-1 do begin cur=where( strcmp(class,cLabel[kClass]), nbCur ) cgPlot, xKey, yKey+kClass*yShi, psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass], /OverPlot cgText, xKey+xShi, yKey+kClass*yShi+yTxt, cLabel[kClass]+'-types: ', color=cColor[kClass] cgText, xKey+xShi*13, yKey+kClass*yShi+yTxt, strtrim(string(nbCur),2), align=1, color=cColor[kClass] endfor cgPS_close, /png, Delete_PS=0 endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- IV --- Classify MOVIS Objects using BDM reference ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--IV.1--Loop Over SSO in MOVIS --------------------------------------------------------- ssoClass = replicate('U',nbSSO) ssoComb = intarr(nbSSO) for kSSO=0, nbSSO-1 do begin distArr=fltarr(nbClass) chi2Arr=fltarr(nbClass) ;--IV.2--Loop Over Taxonomic Classes -------------------------------------------------- for kClass=0, nbClass-1 do begin nbUse = 0 chi2 = 0 ;--IV.3--Loop Over Combination of Filters ------------------------------------------- for kComb1=0, nbComb-1 do begin if ssoCol[kSSO,kComb1] ne -99.99 then $ for kComb2=nbComb-1, kComb1+1, -1 do begin if ssoCol[kSSO,kComb1] ne -99.99 then begin dX = ssoCol[kSSO,kComb1]-covArr[kClass,kComb1,kComb2].meanX dY = ssoCol[kSSO,kComb2]-covArr[kClass,kComb1,kComb2].meanY rotX = dX*covArr[kClass,kComb1,kComb2].cos + dY*covArr[kClass,kComb1,kComb2].sin rotY =-dX*covArr[kClass,kComb1,kComb2].sin + dY*covArr[kClass,kComb1,kComb2].cos chi2 += (rotX/covArr[kClass,kComb1,kComb2].stdA)^2. + $ (rotY/covArr[kClass,kComb1,kComb2].stdB)^2. nbUse++ endif endfor endfor chi2Arr[kClass] = chi2 endfor ;--IV.4--Find Most-Likely Class ------------------------------------------------------- if nbUSe ge 1 then begin order = sort( chi2Arr ) order = order[indgen(6)] ; used = where( ssoCol[kSSO,*] ne -99.99, nbUse ) ; ssoComb[kSSO] = nbUse ; ssoClass[kSSO] = cLabel[order[0]] dof = nbUse-1 print, '' print, kSSO, cLabel[order] ; print, kSSO, ' ['+class[kSSO]+'] ',cLabel[order] print, chi2Arr[order] print, chi2Arr[order]/dof print, dof, 2./dof proba = chisqr_pdf( chi2Arr[order], dof ) print, 1-proba sel = where( (1-proba) ge 0.001, nbSel ) if nbSel gt 0 then begin ssoComb[kSSO] = nbUse ssoClass[kSSO] = cLabel[order[0]] endif endif endfor ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- V --- Selection of Valid Classification ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--V.1--Only Select SSO with More than numSel colors ------------------------------------ sel = where( ssoComb ge numSel, nbSel ) print, 'Number of asteroids with more than '+strtrim(string(numSel),2)+' colors:', nbSel if checkTaxo eq 0 then begin index = where( ssoComb ge numSel and $ ssoCol[*,0] ne -99.99 and $ ssoCol[*,1] ne -99.99 and $ ssoCol[*,2] ne -99.99 and $ ssoCol[*,3] ne -99.99 and $ ssoCol[*,4] ne -99.99 and $ ssoCol[*,5] ne -99.99 ) endif else begin index = indgen(nbSSO) endelse ;--V.2--Validation Figures -------------------------------------------------------------- if showGraphMOVIS eq 1 then begin ;--V.2.1--Loop over Filter Combinations for kComb1=0, nbComb-1 do begin for kComb2=kComb1+1, nbComb-1 do begin ;--V.2.2--Prepare Plot OutLay cgPlot, ssoCol[index,kComb1], ssoCol[index,kComb2], /NoData, xTitle=label[kComb1], yTitle=label[kComb2] ;--V.2.3--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(cLabel[kClass],ssoClass) ) cgPlot, /OverPlot, ssoCol[cur,kComb1], ssoCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--V.2.4--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin cgPlot, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*cos(theta), $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*sin(theta), color=cColor[kClass], /over cgErrPlot, covArr[kClass,kComb1,kComb2].meanX, color=cColor[kClass], $ covArr[kClass,kComb1,kComb2].meanY-covArr[kClass,kComb1,kComb2].stdB, $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB cgErrPlot, covArr[kClass,kComb1,kComb2].meanY, color=cColor[kClass], $ covArr[kClass,kComb1,kComb2].meanX-covArr[kClass,kComb1,kComb2].stdA, $ covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA, /horizontal cgText, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*0.3, $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*0.3, cLabel[kClass], color=cColor[kClass] endfor stop endfor endfor endif ;--V.3--Validation Figures in a File -------------------------------------------------- if createGraphMOVIS eq 1 then begin ;--V.3.1-- Prepare plot if checkTaxo eq 0 then begin filename = 'MOVIS-Result-'+strtrim(string(numSel),2)+'-'+strtrim(string(medCut,format='(F4.2)'),2)+'.eps' endif else begin filename = 'DeMeo-Valid.eps' endelse cgPS_open, Filename=dirTaxo+filename, /metric, /decomposed, /encapsulated, $ xSize=40, ySize=40, language_level=2, /quiet xLen = 0.176 yLen = xLen xSep = 0.01 ySep = 0.01 xBot = 0.085 yBot = 0.05 ;--V.3.2-- Loop over Combinations for kComb1=0, nbComb-1 do begin for kComb2=nbComb-1, kComb1+1, -1 do begin ;--V.3.3--Prepare Plot OutLay ;--V.3.3/A--XY Range lonX = [ covArr[*,kComb1,kComb2].meanX+covArr[*,kComb1,kComb2].stdA, $ covArr[*,kComb1,kComb2].meanX-covArr[*,kComb1,kComb2].stdA, ssoCol[index,kComb1]] lonY = [ covArr[*,kComb1,kComb2].meanY+covArr[*,kComb1,kComb2].stdB, $ covArr[*,kComb1,kComb2].meanY-covArr[*,kComb1,kComb2].stdB, ssoCol[index,kComb2]] ; xRang = [min(lonX), max(lonX)] ; yRang = [min(lonY), max(lonY)] xRang = ccRang[*,kComb1] yRang = ccRang[*,kComb2] ;--V.3.3/B--Graph position pos= [xBot + kComb1*xLen, $ yBot + (kComb2-1)*xLen, $ xBot + (kComb1+1)*xLen, $ yBot + (kComb2 )*xLen ] ;--V.3.3/C--XY Axes On the First column if kComb1 eq 0 then begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTitle=label[kComb2] ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTitle=label[kComb2] endelse ;--V.3.3/D--XY Axes On the other columns endif else begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20) ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20) endelse endelse ;--V.3.4--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(cLabel[kClass],ssoClass) ) cgPlot, /OverPlot, ssoCol[cur,kComb1], ssoCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--V.3.5--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin ;--V.3.5/A--Long-Axis Line xA = [-1,1]*covArr[kClass,kComb1,kComb2].stdA yA = [0,0] rotXA = xA*covArr[kClass,kComb1,kComb2].cos - yA*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYA = xA*covArr[kClass,kComb1,kComb2].sin + yA*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXA, rotYA, color=cColor[kClass], /OverPlot ;--V.3.5/B--Short-Axis Line xB = [0,0] yB = [-1,1]*covArr[kClass,kComb1,kComb2].stdB rotXB = xB*covArr[kClass,kComb1,kComb2].cos - yB*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYB = xB*covArr[kClass,kComb1,kComb2].sin + yB*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXB, rotYB, color=cColor[kClass], /OverPlot ;--V.3.5/C--Ellipse x = covArr[kClass,kComb1,kComb2].stdA * cos(theta) y = covArr[kClass,kComb1,kComb2].stdB * sin(theta) rotX = x*covArr[kClass,kComb1,kComb2].cos - y*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotY = x*covArr[kClass,kComb1,kComb2].sin + y*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotX, rotY, color=cColor[kClass], /OverPlot ; cgPlot, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*cos(theta), $ ; covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*sin(theta), color=cColor[kClass], /over ; ; cgErrPlot, covArr[kClass,kComb1,kComb2].meanX, color=cColor[kClass], $ ; covArr[kClass,kComb1,kComb2].meanY-covArr[kClass,kComb1,kComb2].stdB, $ ; covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB ; cgErrPlot, covArr[kClass,kComb1,kComb2].meanY, color=cColor[kClass], $ ; covArr[kClass,kComb1,kComb2].meanX-covArr[kClass,kComb1,kComb2].stdA, $ ; covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA, /horizontal ; ; cgText, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*0.3, $ ; covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*0.3, cLabel[kClass], color=cColor[kClass] endfor endfor endfor ;--III.4.5-- Display Key pos= [xBot + 3*xLen, $ yBot + (0)*yLen, $ xBot + (5)*xLen, $ yBot + (2 )*yLen ] cgPlot, [0,1], [0,1], /NoData, /NoErase, /Normal, position=pos, xStyle=5, yStyle=5 xKey = 0.2 xShi = 0.03 yKey = 0.90 yShi =-0.075 yTxt = -0.015 for kClass=0, nbClass-1 do begin cur=where( strcmp(ssoClass,cLabel[kClass]), nbCur ) cgPlot, xKey, yKey+kClass*yShi, psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass], /OverPlot cgText, xKey+xShi, yKey+kClass*yShi+yTxt, cLabel[kClass]+'-types: ', color=cColor[kClass] cgText, xKey+xShi*13, yKey+kClass*yShi+yTxt, strtrim(string(nbCur),2), align=1, color=cColor[kClass] endfor cgPS_close, /png, Delete_PS=0 bad = where( strcmp(ssoClass ,'U'),nbu ) print, nbu, nbsso, nbSSO-nbU ; forprint, movis.num, movis.name, ssoClass, ssoComb, textout=2, $ ; format='(I6,1x,A-20,1x,A1,1x,I1)', subset=sel endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- VI --- Classification Validation with DeMeo ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; j2xm: if checkTaxo eq 1 then begin ;ben do comments & all comArr = intarr(nbClass,nbClass)+2 comArr[indgen(nbClass),indgen(nbClass)] = 0 ok_A=[6,7,10] & A=0 & comArr[A,ok_A]=1 & comArr[ok_A,A]=1 ok_B=[2,9] & B=1 & comArr[B,ok_B]=1 & comArr[ok_B,B]=1 ok_C=[9] & C=2 & comArr[C,ok_C]=1 & comArr[ok_C,C]=1 ok_D=[8] & D=3 & comArr[D,ok_D]=1 & comArr[ok_D,D]=1 ok_K=[5,7] & K=4 & comArr[K,ok_K]=1 & comArr[ok_K,K]=1 ok_L=[7] & L=5 & comArr[L,ok_L]=1 & comArr[ok_L,L]=1 ok_Q=[7,10] & Q=6 & comArr[Q,ok_Q]=1 & comArr[ok_Q,Q]=1 ok_S=[10] & S=7 & comArr[S,ok_S]=1 & comArr[ok_S,S]=1 ; ok_T=[8] & T=8 & comArr[T,ok_T]=1 & comArr[ok_T,T]=1 ok_X=[2] & X=9 & comArr[X,ok_X]=1 & comArr[ok_X,X]=1 ok_V=[7] & V=10 & comArr[V,ok_V]=1 & comArr[ok_V,V]=1 xmN = intarr(nbClass,nbClass) xmF = fltarr(nbClass,nbClass) xmL = strarr(nbClass,nbClass) xmV = fltarr(nbClass,nbClass) for kClass=0, nbClass-1 do begin cur = where( strcmp( Class, cLabel[kClass]), nbValid ) local = ssoClass[cur] uniqClass = local[ uniq( local, sort(local) ) ] nbUniq=n_elementS(uniqClass) numPerClass=intarr(nbUniq) frqPerClass=fltarr(nbUniq) for kU=0, nbUniq-1 do begin pp = where( strcmp(local,uniqClass[kU]), nbPP ) numPerClass[kU]=nbPP endfor frqPerClass=100.*numPerClass/nbValid ord=reverse(sort(numPerClass)) if strcmp( cLabel[kClass], 'T' ) then begin ord=sort(numPerClass) endif xmN[kClass,0:nbUniq-1] = numPerClass[ord] xmF[kClass,0:nbUniq-1] = frqPerClass[ord] xmL[kClass,0:nbUniq-1] = uniqClass[ord] print, '' print, cLabel[kClass], ' -- ', uniqClass[ord[0]], numPerClass[ord[0]], frqPerClass[ord[0]], $ format='(A1,A4,A1,3x,I3,3x,F6.2)' if nbUniq gt 1 then begin for kU=1, nbUniq-1 do begin print, ' ', ' ', uniqClass[ord[kU]], numPerClass[ord[kU]], frqPerClass[ord[kU]], $ format='(A1,A4,A1,3x,I3,3x,F6.2)' endfor print, ' ', ' ', 'Tot', nbValid, total(frqPerClass), $ format='(A1,A3,A3,2x,I3,3x,F6.2)' endif endfor statValid=fltarr(nbClass,3) statValid[*,0] = xmF[*,0] cgPS_open, Filename=dirTaxo+'DeMeo-Valid-Fraction.eps', /metric, /decomposed, /encapsulated, $ xSize=25, ySize=18, language_level=2, /quiet cgBarPlot, reverse(xmF[*,0]), /Rotate, /normal, $ xStyle=1, xRange=[0,100], $ yStyle=1, yRange=[0,nbClass-1], $ barNames=reverse(cLabel), colors=reverse(cColor), $ xTitle='Fraction of classification', $ yTitle='Euclid classification', $ barCoords=yArr, $ position=[0.08,0.12,0.97,0.97] barSize=0.6675 barOffset=-0.33 last = xmF[*,0] for kE=1, nbClass-1 do begin colorVec=replicate('Black',nbClass) ; for kClass=nbClass-1,0,-1 do begin for kClass=0, nbClass-1 do begin p=where( strcmp(cLabel, xmL[kClass,kE]) ) if p[0] ne -1 then begin colorVec[kClass] = cColor[p[0]] if comArr[kClass,p[0]] lt 2 then xmV[kClass,kE]=1 statValid[*, comArr[kClass,p[0]]] += xmF[*,kClass] if xmV[kClass,kE] eq 1 and xmF[kClass,kE] gt 0.5 then begin base = last last[kClass] += xmF[kClass,kE] xx = base[kClass] + [0,xmF[kClass,kE],xmF[kClass,kE],0,0] yy = [0,0,1,1,0]*barSize + barOffSet + yArr[nbClass-1-kClass] cgPlot, /over, xx, yy, color=colorVec[kClass] cgText, mean(xx[0:1]), yy[1] + (yy[2]-yy[1])/4., align=0.5, xmL[kClass,kE], color=colorVec[kClass] endif endif endfor endfor cgPS_close, /png, Delete_PS=0 endif end