;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-What to do genSUM = 0 genPLOT= 01 ;---------- 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/' dirLC = dirEuclid+'lc/' ;--- survey limits vLimit=24. elongation = 90. ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Generate a H/diameter vs period file -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if genSUM eq 1 then begin ;--II.1-- Read PDS LC summary lc = pds_lcparam_getcatalog('*') ao = readAstOrb() nbLC = n_elements(lc) h = fltarr(nbLC) ;--II.2-- Numbered and non-numbered objects wn = where( lc.num ne 0, nwn, comp=wo, ncomp=nwo ) h[wn] = ao[wn].h for k=0, nwo-1 do begin p=where( strcmp(lc[wo[k]].des, ao.name, /fold)) if p[0] ne -1 then h[wo[k]] = ao[p[0]].h endfor ;--II.3-- Clean Sample from per=0 and h=0 g = where( h ne 0 and lc.p.val gt 0, nbSSO) h=h[g] lc=lc[g] ;--II.4-- Find dynamical class dyn=strarr(nbSSO) wn=where( lc.num ne 0, nwn, comp=wo, ncomp=nwo ) for k=0, nwn-1 do begin class = astClass(lc[wn[k]].num) dyn[wn[k]]= class.type endfor for k=0, nwo-1 do begin class = astClass(lc[wo[k]].des) dyn[wo[k]]= class.type endfor ;-II.5-- Export results on disk good = where( ~strcmp(lc.p.flag, 'u',/fold) ) forprint, lc.num, lc.des, dyn, h, $ lc.p.val, lc.p.qual, lc.p.flag, lc.p.note, $ lc.a.min, lc.a.max, lc.a.flag, $ format='(F7,1x,A-17,1x,A-7,1x,F5.2,1x,'+$ 'F14.8,1x,A2,1x,A1,1x,A-15,1x,'+$ 'F4.2,1x,F4.2,1x,A1)', $ ; subset=good, $ textout=dirLC+'lc4euclid.dat', /NoComment forprint, lc.num, lc.des, dyn,h, $ lc.p.val, lc.p.qual, lc.p.flag, lc.p.note, $ lc.a.min, lc.a.max, lc.a.flag, $ format='(I7," , ",A-17," , ",A-12," , ",F5.2," , ",'+$ 'F14.8," , ",A-2," , ",A1," , ",A-15," , ",'+$ 'F6.2," , ",F6.2," , ",A1)', $ comment='Num, Des, Dyn, H, P, P qual, P flag, P note, A min, A max, A flag', $ ; subset=good, $ textout=dirLC+'lc4euclid.csv' spawn, "awk 'NR>1' "+dirLC+'comet.csv >> '+dirLC+'lc4euclid.csv' endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Generate the H vs P plot -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if genPLOT eq 1 then begin ;--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, popS, /Silent popName = strtrim(popName,2) nbPop=n_elements(popName) sel = where( popMain eq 1, nbSel ) popName = popName[sel] popC = popC[sel] popS = popS[sel] nbPop = nbSel print, 'Populations considered: ('+strtrim(string(nbPop),2)+'): ', popName ;--III.2-- Read the H/D/P file readcol, dirLC+'lc4euclid.csv', num, des, dyn, h, $ pval, pqual, pflag, pnote, amin, amax, aflag, delimiter=',', $ format='(I,A,A,A,F,A,A,A,F,F,A)', /silent pqual=strtrim(pqual,2) dyn=strtrim(dyn,2) hun=where( strCmp(dyn,'Hungari') ) dyn[hun]='MB' vecQual = intarr(n_elements(num)) selQ = where( strCmp(pqual,'2',1) or strCmp(pqual,'3',1) ) vecQual[selQ] = 1 ;--III.3-- Plot parameters and open xRang = [33,2] yRang = [1e4,1e-2] xShi = -0.5 xPos = 31.7 + indgen(nbPop)*xShi kPos=0 xKey = 8 yKey = cgLogGen( nbPop, start=7e-1, finish=2.5e-2) yKey2 = cgLogGen( nbPop, start=8e-1, finish=3e-2) ; yKey =yKey [[2,0,1,3,4,5]] ; yKey2=yKey2[[2,0,1,3,4,5]] cgPS_open, Filename=dirLC+'Euclid-LC.eps', /metric, /decomposed, /encapsulated, $ xSize=35, ySize=18.9, language_level=2, /quiet ;--III.4-- Plot xTickName = strtrim(string( 32 - indgen(20)*2, format='(I2)'),2) xTickName[0:1]=' ' cgPlot, h ,pval, /yLog, /Normal, /NoData, $ xStyle=1+8, xRange=xRang, xTitle='Absolute magnitude', $ yStyle=1, yRange=yRang, yTitle='Rotation period (h)', $ xTickInterval=2, xMinor=2, xTickName=xTickName, $ yTickFormat='exponent', $ psym='Open circle', $ symSize=0.5, color='Light gray', $ position=[0.08,0.10,0.95,0.95] dRang = (1329/sqrt(0.20)) * 10^(-0.2 * xRang) cgAxis, xAxis=1, xRange=dRang, /xLog, xStyle=1, xMinor=10, $ xTickName=reverse([' ','1000','100','10','1','0.1',' ']) cgText, /Normal, 0.85,0.967, 'km' cgText, /Normal, 0.10,0.967, 'Diameter' ; cgAxis, 28, yAxis=1, yRange=yRang, /yLog, yStyle=1, yTickName=replicate(' ',10) cgAxis, 28, yAxis=0, yRange=yRang, /yLog, yStyle=1, yTickName=replicate(' ',10) ; print, dyn[ uniq(dyn, sort(dyn)) ] stat=fltarr(nbPop,5) mu=fltarr(nbPop) sig=fltarr(nbPop) ; for kPop=0, nbPop-1 do begin ord=[2,0,1,3,4,5,6] for kOrd=0, nbPop-1 do begin kPop = ord[kOrd] cur = where( dyn eq popName[kPop] and vecQual eq 1 and h lt 28, nbCur ) ; print, popName[kPop], nbCur if cur[0] ne -1 then begin if strcmp(popName[kPop],'Comet') then pval *= 24. stat[kPop,*] = dataBox(pval[cur]) mu[kPop] = mean(pval[cur]) sig[kPop] = stddev(pval[cur]) print, popName[kPop], reform(stat[kPop,*]), nbCur, $ format='(A-10,2x,F9.3,F9.3,F9.3,F9.3,F9.3,2x,I5)' ; cgPlot, /OverPlot, h[cur] ,pval[cur], color=popC[kPop], psym='Filled circle' ; cgPlot, /OverPlot, xKey+0.3 ,yKey[kPop], color=popC[kPop], psym='Filled circle' cgPlot, /OverPlot, h[cur] ,pval[cur], color=popC[kPop], psym=popS[kPop] cgPlot, /OverPlot, xKey+0.3 ,yKey[kPop], color=popC[kPop], psym=popS[kPop] cgText, xKey ,yKey2[kPop], popName[kPop] cgText, xKey-5.0 ,yKey2[kPop], strtrim(string(nbCur,format='(I)'),2), align=1 cgErrPlot, xPos[kPos], stat[kPop,0],stat[kPop,4], color=popC[kPop] bLen=0.25 boxX = bLen*0.5*[-1,1,1,-1,-1] + xPos[kPos] boxY = [ stat[kPop,1], stat[kPop,1], stat[kPop,3], stat[kPop,3], stat[kPop,1] ] cgColorFill, boxX, boxY, color=popC[kPop] cgPlot, /OverPlot, boxX, boxY cgPlot, /OverPlot, bLen*0.5*[-1,1]+xPos[kPos], [ stat[kPop,2], stat[kPop,2] ] kPos++ endif endfor cgPS_close, /PNG, delete_PS=0 endif end