;goto, j2Plot ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-What to do doCOMP = 0 doFIG = 1 ;---------- Directories spawn, 'hostname', host case strTrim(host) of 'hyperion': dirEuclid = '/data/euclid/' 'endymion': dirEuclid = '/home/bcarry/work/data/euclid/' else: stop endcase dirFilter = dirEuclid + 'filters/' dirTaxo = dirEuclid + 'taxonomy/' dirStat = dirEuclid + 'stat/' dirSSO = dirEuclid + 'ssos/astermpi/' dirRate = dirEuclid + 'rate/' dirExten = dirEuclid + 'exten/' dirRS = dirEuclid + 'rs/' ;--- survey limits maxEc = 30 limEcl = 15 limGal = 20. limElong = 11.5 band = 5. mjd2000 = date_conv( '2000-01-01T00:00:00', 'JULIAN' ) label=['VIS','NISP','Y','J','H'] DIT = [565, 565, 121, 116, 81.] ;-s pix = [0.1, 0.3, 0.3, 0.3, 0.3] ;-arcsec ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- SSO populations, time, and RS -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Read SSO populations readcol, dirEuclid+'pop.rate', 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 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Compute Apparent Rate Statistics -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if doCOMP eq 1 then begin xRang=[0,100] aRang=[-90.,90] binRate = 0.1 ; "/h binAngle= 5. ; degree nbBinRate = (xRang[1]-xRang[0])/binRate + 1 histRate = intarr(nbPop, nbBinRate) statRate = fltarr(nbPop,5) stdRate = fltarr(nbPop,2) nbBinAngle = (aRang[1]-aRang[0])/binAngle + 1 histAngle = intarr(nbPop, nbBinAngle) statAngle = fltarr(nbPop,5) for kPop=0, nbPop-1 do begin ;--V.2.1-- Define file suffix for each population case popName[kPop] of 'NEA': suffix ='nea' 'MC': suffix = 'mars_crosser' 'MB': suffix = 'mb' 'Trojan': suffix = 'trojan' 'Centaur': suffix = 'centaur' 'KBO': suffix = 'kbo' 'Comet': suffix = 'comet' endcase nameEPH = 'ephem_'+suffix+'.csv' readcol, dirSSO+nameEPH, delimiter=',', /Silent, $ id, JD, RA, DEC, Dobs, MagV, Phase, SunElong, Dh, dRA, dDEC, dDobs, Class, Name, $ format='(L,D,F,F,F,F,F,F,F,F,F,F,A,A)' ; sel = where( sunElong ge 87 and sunElong le 110, nbSel) sel = where( sunElong ge 88.5 and sunElong le 91.5, nbSel) mu = sqrt( (dRA[sel]*cos(dec[sel]*!DTOR))^2. + dDec[sel]^2. ) mu *= (3600 / 24.) histoMu = histogram( mu, binSize=binRate, min=xRang[0], max=xRang[1], loc=rArr ) histRate[kPop,*] = histoMu statRate[kPop,*] = dataBox(mu) stdRate[kPop,*] = [mean(mu), stddev(mu)] print, popName[kPop], nbSel/15., statRate[kPop,*], format='(A-10,2x,I6,2x,5(F7.2,2x))' endfor rate = statRate[*,2] ; [0.1,1,10, 50, 100.] nRate = n_elements( rate) openw, idEx, dirRate+'euclid-rate.csv',/get_lun for kRate=0, nRate-1 do $ printf, idEx, popName[kRate], $ statRate[kRate,0], statRate[kRate,1], statRate[kRate,2], statRate[kRate,3], statRate[kRate,4], $ statRate[kRate,2] * ( DIT/pix/3600.), $ format='(A-10,", ",F6.1,", ",F6.1,", ",F6.1,", ",F6.1,", ",F6.1,'+$ '", ",F6.1," , ",F6.1," , ",F6.1," , ",F6.1," , ",F6.1)' free_lun, idEx openw, idEx, dirRate+'euclid-rate.tex',/get_lun for kRate=0, nRate-1 do $ printf, idEx, popName[kRate], $ statRate[kRate,2], statRate[kRate,1]-statRate[kRate,2], statRate[kRate,3]-statRate[kRate,2], $ statRate[kRate,2] * ( DIT/pix/3600.), $ format='(A-10," & ",F4.1,"$_{\scriptscriptstyle ",F5.1,"}^{\scriptscriptstyle ",F+5.1,'+$ '"}$ & ",F5.1," & ",F5.1," & ",F5.1," & ",F5.1," & ",F5.1," \\")' free_lun, idEx endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Build Graphics -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if doFIG eq 1 then begin if doCOMP ne 1 then begin readcol, dirRate+'euclid-rate.csv', delimiter=',', /Silent, $ popi, sr0, sr1, sr2, sr3, sr4, $ pix1, pix2, pix3, pix4, pix5, $ format='(A,F,F,F,F,F,F,F,F,F,F)' statRate = fltarr(nbPop,5) statRate[*,0] = sr0 statRate[*,1] = sr1 statRate[*,2] = sr2 statRate[*,3] = sr3 statRate[*,4] = sr4 tooLow = where( statRate le 0.1 ) statRate[tooLow]=0.1 tooHigh = where( statRate ge 1000 ) statRate[tooHigh]=1000 endif cgPS_open, Filename=dirRate+'Euclid-Rate.eps', /metric, /decomposed, /encapsulated, $ xSize=35, ySize=12, language_level=2, /quiet cgPlot, 0,0, /NoErase, /NoData, /Normal, $ xStyle=1, yRange=[0,nbPop+1], $ yStyle=1, xRange=[0.1,1000], /xLog, $ position=[0.10,0.18,0.96,0.98], $ xTickInterval=30, xMinor=10, $ yTickInterval=1, yMinor=1, $ xTitle='Apparent rate ("/h)', $ yTickName=[' ',popName,' '], $ xTickName=['0.1','1','10','100','1000'], $ xTickLen=0.05 for kPop=0, nbPop-1 do begin errLow = statRate[kPop,2]-statRate[kPop,0] errHigh= statRate[kPop,4]-statRate[kPop,2] if errLow lt 0.1 then errLow = 0.1 if errHigh gt 1000 then errHigh = 1000 ; cgErrPlot, kPop+1, statRate[kPop,0], statRate[kPop,4], color=popC[kPop], /Horizontal, $ ; err_width=0.2 cgPlot, /OverPlot, statRate[kPop,2], kPop+1, $ err_xLow=errLow, $ err_xHigh=errHigh,$ color=popC[kPop], err_width=0.025 print, reform(statRate[kPop,*]) bLen=0.35 boxX = [ statRate[kPop,1], statRate[kPop,1], statRate[kPop,3], statRate[kPop,3], statRate[kPop,1] ] boxY = bLen*0.5*[-1,1,1,-1,-1] + kPop+1 cgColorFill, boxX, boxY, color=popC[kPop] cgPlot, /OverPlot, boxX, boxY if strCmp(popName[kPop],'MB') then col='white' else col='black' cgPlot, /OverPlot, [ statRate[kPop,2], statRate[kPop,2] ], bLen*0.5*[-1,1]+kPop+1, color=col endfor cgPS_close, /PNG, Delete_PS=0 endif end