;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;----What to do ;---------- Directories spawn, 'hostname', host case strTrim(host) of 'hyperion': dirEuclid = '/observ/euclid/' 'endymion': dirEuclid = '/home/bcarry/work/data/euclid/' else: stop endcase dirFilter = dirEuclid+'filters/' dirTaxo = dirEuclid+'taxonomy/' dirStat = dirEuclid + 'stat/' dirFrac = dirEuclid + 'fraction/' dirPhase= dirEuclid + 'phase/' ;--- survey limits limEcl = 15 limGal = 20. vLimit = 24.5 ;-- YJH -- 24.5 in V, but smearing... binH = 0.25 ;-Geometry of observation euclid = {min:87., mean:90., max:110., shift:+0.0, symb:'Filled circle', size:1} lsst = {min:45., mean:90., max:180., shift:-0.1, symb:'Open square', size:0.75} gaia = {min:45., mean:90., max:145., shift:+0.1, symb:'Open diamond', size:1.0} ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Asteroid Phase Curves -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-General phase phase={step:0.5, min:0., max:120., N:0, arr:ptr_new(/ALLOCATE_HEAP)} phase.N = (phase.max-phase.min)/phase.step + 1 *phase.arr = phase.min + phase.step*findgen(phase.N) ;--- Asteroid phase curve ESC = replicate({num:0, class:'', A:0.,D:0., k:0., f:ptr_new(/ALLOCATE_HEAP),style:0, col:''},3) ESC[0] = {num:317,class:'E',A:0.26,D:0.08,k:-0.47,f:ptr_new(/ALL), style:0, col:'Cornflower blue'} ;-317 ESC[1] = {num:270,class:'S',A:0.50,D:0.09,k:-0.62,f:ptr_new(/ALL), style:1, col:'Red'} ;-270 ESC[2] = {num:313,class:'C',A:1.00,D:0.18,k:-0.75,f:ptr_new(/ALL), style:2, col:'Gray'} ;-747 ESC[*].D /= !DTOR ESC[*].k *= !DTOR *esc[0].f = ssoPhase(ESC[0],*phase.arr,/LCI) *esc[1].f = ssoPhase(ESC[1],*phase.arr,/LCI) *esc[2].f = ssoPhase(ESC[2],*phase.arr,/LCI) *esc[0].f /= max( *esc[0].f ) *esc[1].f /= max( *esc[1].f ) *esc[2].f /= max( *esc[2].f ) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- SSO populations -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Read SSO populations readcol, dirEuclid+'pop', delimiter=',', $ format='(A,I1,F4.1,F4.1,F4.1,F4.1,F3.1,F3,A,A)', $ popName, popMain, popHmin, popHmax, xMin, xMax, yMin, yMax, popC, poS, /Silent popName = strtrim(popName,2) nbPop=n_elements(popName) popPrint = where( popMain eq 1, nbPopPrint ) print, 'Populations considered: ('+strtrim(string(nbPop),2)+'): ', popName group1 = ['Comet', 'NEA','MC', 'MB'] group2 = ['Trojan', 'Centaur','KBO'] nbG1 = n_elements(group1) nbG2 = n_elements(group2) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Make Plots -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Plot definition -------------------------------------------------------------------- xBot = 2500 yBot = 2000 yLen1 = 6000 yLen2 = 15000 ySep = 1200 xLen1 = 19000. xLen2 = 13000 xLen3 = 6000 xSep = 150 yShift = 1 pos1 = [ xBot, yBot, xBot+xLen1, yBot+yLen1 ] pos2 = [ xBot, yBot+yLen1+ySep, xBot+xLen2, yBot+yLen1+ySep+yLen2 ] pos3 = [ xBot+xLen2+xSep, yBot+yLen1+ySep, xBot+xLen1, yBot+yLen1+ySep+yLen2 ] cgPS_open, Filename=dirPhase+'Euclid-Phase.eps', /metric, /decomposed, /encapsulated, $ xSize=22, ySize=25.5, language_level=2, /quiet ;--IV.2-- Different classes ------------------------------------------------------------------ cgPlot, 0,0, /NoData, /NoErase, /Device, $ xStyle=1, xRange=[0,90], xTickLen=0.05, $ yStyle=1, yRange=[0,1], yTickLen=0.02, $ yTickInt=0.2, yMinor=2, $ xTitle='Phase angle (!Uo!N)', $ yTitle='Relative flux', $ position=pos1 xKey = 60 yKey = 0.8 yKS = -0.125 xKL = 10 for k=0,2 do begin cgPlot, /OverPlot, *phase.arr, *esc[k].f, color=esc[k].col, linestyle=esc[k].style cgPlot, /OverPlot, xKey+[0,xKL], [yKey,yKey]+k*yKS, color=esc[k].col, linestyle=esc[k].style cgText, xKey+xKL*1.2, yKey+k*yKS-0.05, esc[k].class+'-type' endfor ;--IV.3-- Plot definition -------------------------------------------------------------------- cgPlot, 0,0, /NoData, /NoErase, /Device, $ xStyle=1+8, xRange=[0,90], $ yStyle=1+8, yRange=[0,4], $ xTickInt=20, xMinor=4, $ yTickInt=1 , yMinor=5, yTickLen=0.03, yTickFormat='(F3.1)', $ yTitle='Relative flux (shifted for clarity)', $ position=pos2 xKey = 5 ySurv = 0.1 modP = 4 for k=0, nbG1-1 do begin kPop=where( popName eq group1[k] ) print, popName[kPop] ;--IV.3.1-- general S type phase function cgPlot, /OverPlot, *phase.arr, *esc[1].f + k*yShift, color=popC[kPop] cgText, xKey, max(*esc[1].f) + k*yShift, popName[kPop], color=popC[kPop] ;--IV.3.2-- Euclid Survey range trash = ssoHtoV(pop=popName[kpop], elong=euclid.Min , geometry=gMin ) trash = ssoHtoV(pop=popName[kpop], elong=euclid.Mean, geometry=gAvg ) trash = ssoHtoV(pop=popName[kpop], elong=euclid.Max , geometry=gMax ) allPhase = [gMin.phase, gAvg.phase, gMax.phase] minPhase = min(allPhase,/Nan) maxPhase = max(allPhase,/Nan) if popName[kPop] eq 'Centaur' then maxPhase = 9 range = where( *phase.arr ge minPhase and *phase.arr le maxPhase and (indgen(phase.N) mod modP) eq 1) cgPlot, /OverPlot, (*phase.arr)[range], psym=euclid.symb, symsize=euclid.size, $ (*esc[1].f)[range] + k*yShift, color=popC[kPop] ;--IV.3.3-- LSST Survey range trash = ssoHtoV(pop=popName[kpop], elong=lsst.Min , geometry=gMin ) trash = ssoHtoV(pop=popName[kpop], elong=lsst.Mean, geometry=gAvg ) trash = ssoHtoV(pop=popName[kpop], elong=lsst.Max , geometry=gMax ) allPhase = [gMin.phase, gAvg.phase, gMax.phase] minPhase = min(allPhase,/Nan) maxPhase = max(allPhase,/Nan) if popName[kPop] eq 'Centaur' then maxPhase = 9 range = where( *phase.arr ge minPhase and *phase.arr le maxPhase and (indgen(phase.N) mod modP) eq 1) ; range = where( *phase.arr ge minPhase and *phase.arr le maxPhase ) cgPlot, /OverPlot, (*phase.arr)[range], psym=lsst.symb, symsize=lsst.size, $ (*esc[1].f)[range] + k*yShift-ySurv, color=popC[kPop] ;--IV.3.4-- Gaia Survey range trash = ssoHtoV(pop=popName[kpop], elong=gaia.Min , geometry=gMin ) trash = ssoHtoV(pop=popName[kpop], elong=gaia.Mean, geometry=gAvg ) trash = ssoHtoV(pop=popName[kpop], elong=gaia.Max , geometry=gMax ) allPhase = [gMin.phase, gAvg.phase, gMax.phase] minPhase = min(allPhase,/Nan) maxPhase = max(allPhase,/Nan) if popName[kPop] eq 'Centaur' then maxPhase = 9 range = where( *phase.arr ge minPhase and *phase.arr le maxPhase and (indgen(phase.N) mod modP) eq 1) ; range = where( *phase.arr ge minPhase and *phase.arr le maxPhase ) cgPlot, /OverPlot, (*phase.arr)[range], psym=gaia.symb, symsize=gaia.size, $ (*esc[1].f)[range] + k*yShift+ySurv, color=popC[kPop] endfor ;--IV.1-- Plot definition -------------------------------------------------------------------- cgPlot, 0,0, /NoData, /NoErase, /Device, $ xStyle=1+8, xRange=[0,15], $ yStyle=1+8, yRange=[0,5], $ xTickInt=10, xMinor=5, $ yTickInt=1 , yMinor=5, yTickLen=0.06, $;yTickFormat='(F3.1)', $ yTickName=replicate(' ',10), $ position=pos3 xKey = 5 ySurv = 0.1 for k=0, nbG2-1 do begin kPop=where( popName eq group2[k] ) print, popName[kPop] ;--IV.3.1-- general S type phase function cgPlot, /OverPlot, *phase.arr, *esc[1].f + k*yShift, color=popC[kPop] cgText, xKey, max(*esc[1].f) + k*yShift, popName[kPop], color=popC[kPop] ;--IV.3.2-- Euclid Survey range trash = ssoHtoV(pop=popName[kpop], elong=euclid.Min , geometry=gMin ) trash = ssoHtoV(pop=popName[kpop], elong=euclid.Mean, geometry=gAvg ) trash = ssoHtoV(pop=popName[kpop], elong=euclid.Max , geometry=gMax ) allPhase = [gMin.phase, gAvg.phase, gMax.phase] minPhase = min(allPhase,/Nan) maxPhase = max(allPhase,/Nan) if popName[kPop] eq 'Centaur' then maxPhase = 9 range = where( *phase.arr ge minPhase and *phase.arr le maxPhase ) cgPlot, /OverPlot, (*phase.arr)[range], psym=euclid.symb, symsize=euclid.size, $ (*esc[1].f)[range] + k*yShift, color=popC[kPop] ;--IV.3.3-- LSST Survey range trash = ssoHtoV(pop=popName[kpop], elong=lsst.Min , geometry=gMin ) trash = ssoHtoV(pop=popName[kpop], elong=lsst.Mean, geometry=gAvg ) trash = ssoHtoV(pop=popName[kpop], elong=lsst.Max , geometry=gMax ) allPhase = [gMin.phase, gAvg.phase, gMax.phase] minPhase = min(allPhase,/Nan) maxPhase = max(allPhase,/Nan) if popName[kPop] eq 'Centaur' then maxPhase = 9 range = where( *phase.arr ge minPhase and *phase.arr le maxPhase ) cgPlot, /OverPlot, (*phase.arr)[range], psym=lsst.symb, symsize=lsst.size, $ (*esc[1].f)[range] + k*yShift-ySurv, color=popC[kPop] ;--IV.3.4-- Gaia Survey range trash = ssoHtoV(pop=popName[kpop], elong=gaia.Min , geometry=gMin ) trash = ssoHtoV(pop=popName[kpop], elong=gaia.Mean, geometry=gAvg ) trash = ssoHtoV(pop=popName[kpop], elong=gaia.Max , geometry=gMax ) allPhase = [gMin.phase, gAvg.phase, gMax.phase] minPhase = min(allPhase,/Nan) maxPhase = max(allPhase,/Nan) if popName[kPop] eq 'Centaur' then maxPhase = 9 range = where( *phase.arr ge minPhase and *phase.arr le maxPhase ) cgPlot, /OverPlot, (*phase.arr)[range], psym=gaia.symb, symsize=gaia.size, $ (*esc[1].f)[range] + k*yShift+ySurv, color=popC[kPop] endfor cgPS_close, /PNG, delete_PS=0 stop ; cgPlot, phase, f0, $ ; xTitle='Phase angle (!Do!N)', $ ; yTitle='Relative flux' ; cgPlot, /OverPlot, phase, f1, color='red', linestyle=2 ; cgPlot, /OverPlot, phase, f2, color='blue' ;;stop ;; m0 = 2.5*alog10(f0) ;; m1 = 2.5*alog10(f1) ;; m2 = 2.5*alog10(f2) ;; ;;; window, 2 ;; cgPlot, phase, m0-max(m0), yRange = [-6,0.5] ;; cgPlot, /OverPlot, phase, m1-max(m1), color='red', linestyle=2 ;; cgPlot, /OverPlot, phase, m2-max(m2), color='blue' ; ; ; ; ; cgPlot, phase, f1/max(f1), $ ; xTitle='Phase angle (!Do!N)', $ ; yTitle='Relative flux', $ ; xRange=[0.5, 60], $ ; yRange=[0, 7] yShift = 1 ;--IV.2-- Loop over Populations kPlot=0 for kPop=0, nbPop-1 do begin if popMain[kPop] eq 1 then begin ; trash = ssoHtoV(pop=popName[kpop], elong=euclid.Min , geometry=gMin ) ; trash = ssoHtoV(pop=popName[kpop], elong=euclid.Mean, geometry=gAvg ) ; trash = ssoHtoV(pop=popName[kpop], elong=euclid.Max , geometry=gMax ) trash = ssoHtoV(pop=popName[kpop], elong=lsst.Min , geometry=gMin ) trash = ssoHtoV(pop=popName[kpop], elong=lsst.Mean, geometry=gAvg ) trash = ssoHtoV(pop=popName[kpop], elong=lsst.Max , geometry=gMax ) allPhase = [gMin.phase, gAvg.phase, gMax.phase] minPhase = min(allPhase,/Nan) maxPhase = max(allPhase,/Nan) if popName[kPop] eq 'Centaur' then maxPhase = 9 ; nbPhase = (maxPhase-minPhase)/newStep + 1 ; newPhase = minPhase + newStep * findgen(nbPhase) ; ; newF1 = ssoPhase(ESC[1],newPhase) ; maxNew = max( newF1, pMax ) ; maxRef = f1( where( phase eq newPhase[pMax] ) ) ; print, phase( where( phase eq newPhase[pMax] ) ) ; print, newPhase[pMax] ; newF1 *= maxRef/maxNew ; cgPlot, /OverPlot, phase, f1+kPlot*yShift ; cgPlot, /OverPlot, newPhase, newF1 +kPlot*yShift, psym='Open circle' ; print, popName[kpop], minPhase, maxPhase, format='(A-20,2x,I3,2x,I3)' ;; cgText, 50, ssoPhase(ESC[1],50)+kPlot*yShift, popName[kpop] kPlot++ endif endfor end