;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; readData = 01 figOnly = 0 binSize = 1 kMain=0 if figOnly eq 1 then goto, j2fig ;---------- 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/' dirRate = dirEuclid + 'rate/' dirFrac = dirEuclid + 'fraction/' dirSSO = dirEuclid + 'ssos/' dirExten = dirEuclid + 'exten/' dirRS = dirEuclid + 'rs/' ;--- survey limits maxEc = 30 limEcl = 15 limGal = 20. band = 5. fov =0.5 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- SSO populations -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if readData eq 1 then begin ;--II.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) popPrint = where( popMain eq 1, nbPopPrint ) print, 'Populations considered: ('+strtrim(string(nbPop),2)+'): ', popName ;--II.2-- Read expected number of SSO at Euclid limiting mag (from euclid-stat.pro) readcol, dirStat+'euclid-stat.csv', delimiter=',', /Silent, $ statName, statF, statMin, stat25, stat50, stat75, statMax, format='(A,F,F,F,F,F,F)' ; totalNumber = total(nbAvail) ;--II.3-- Read expected number of SSO at Euclid limiting mag (from euclid-stat.pro) readcol, dirStat+'euclid-disco.csv', delimiter=',', /Silent, $ discoName, discoF, discoMin, disco25, disco50, disco75, discoMax, format='(A,F,F,F,F,F,F)' ;--II.4-- Read ASTORB with classes through Discovery Table ; ssoNow = ssoDiscovery_read() ;--II.5-- Read Euclid times readcol, dirEuclid+'time.sso', ep, format='(A)', /silent nbEp = n_elements( ep ) kEp=0 endif ; nbPop=5 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Calibration Fields -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Read Calibration file fileCal = 'calib.csv' cal = euclidReadRS( dirRS+fileCal ) nbCal = n_elements( cal.id ) ;--III.2-- Distribution of CalField over Ecliptic Latitude cgHistoPlot, cal.lat, minInput=0, maxInput=90, binSize=binSize, loc=iArr, histData=histoCal ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Compute the number of SSO toward Ecliptic plane -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ; cgHistoPlot, ssoNow.i, minInput=0, maxInput=90, binSize=1., loc=iArr, histData=histoI ssoHisto = fltarr( nbPop, n_elements(iArr) ) calHisto = fltarr( nbPop, n_elements(iArr) ) obsCal = fltarr( nbPop, n_elements(iArr) ) for kPop=0, nbPop-1 do begin ; kPop=5 ;--III.1-- Define file suffix for each population case popName[kPop] of 'NEA': suffix ='nea' 'Aten': suffix = 'nea-aten' 'Apollo': suffix = 'nea-apollo' 'Amor': suffix = 'nea-amor' 'MC': suffix = 'mcs' 'MB': suffix = 'mba' 'Hungaria': suffix = 'mba-hun' 'IMB': suffix = 'mba-imb' 'MMB': suffix = 'mba-mmb' 'OMB': suffix = 'mba-omb' 'Cybele': suffix = 'mba-cyb' 'Hilda': suffix = 'mba-hil' 'Trojan': suffix = 'jta' 'Centaur': suffix = 'cen' 'KBO': suffix = 'kbo' 'Classical': suffix = 'kbo-classical' 'SDO': suffix = 'kbo-sdo' 'Comet': suffix = 'com' endcase ;--III.2-- Build histogram of inclination if popMain[kPop] eq 1 then begin ;--III.2.1-- Read SSO RA/Dec long/lat at current epoch if suffix ne 'com' then $ readcol, dirSSO+suffix+'-'+ep[kEp]+'.csv', delimiter=',', /Silent, $ id1, ra1, dec1, d1, dra1, ddec1, h1, g1, class1, name1, glon1, glat1, elong1, elat1, $ format='(I,F,F,F,F,F,F,F,A,A,F,F,F,F)' $ else $ readcol, dirSSO+suffix+'-'+ep[kEp]+'.csv', delimiter=',', /Silent, $ id1, ra1, dec1, d1, dra1, ddec1, h1, r1, d1, a1,a2,a3,class1, iau1, name1, glon1, glat1, elong1, elat1, $ format='(I,F,F,F,F,F,F,F,F,F,F,F,A,A,A,F,F,F,F)', skipline=1 cgHistoPlot, elat1, minInput=0, maxInput=90, binSize=binSize, loc=iArr, histData=histoI ssoHisto[kPop,*] = float(histoI)/total(histoI) calHisto[kPop,*] = 0.5 * histoCal * reform(ssoHisto[kPop,*]) / (360. * cos(iArr*!DTOR) ) ;stop endif endfor wdelete ;stop j2fig: ;--III.3-- Plot SSO Distribution over Latitude cgPS_open, Filename=dirExten+'Euclid-Latitude.eps', /metric, /decomposed, /encapsulated, $ xSize=35, ySize=24, language_level=2, /quiet xRang=[0,70] cgPlot, 0,0, /NoData, /Normal, $ xStyle=1, xRange=xRang, $ yStyle=1, yRange=[1, 1e7], /yLog, $ xTitle='Ecliptic latitude (!Uo!N)', $ yTitle='Number of SSOs per latitude degree', $ xTickInterval=5, xMinor=5, $ yTickInterval=1, yMinor=5, $ position=[0.09,0.105,0.96,0.80] for kPop=0, nbPop-1 do begin if popMain[kPop] eq 1 then begin case popName[kPop] of 'Centaur': begin mid = stat75[kPop] low = stat50[kPop] hig = statMax[kPop] end 'Comet': begin mid = stat75[kPop] low = stat50[kPop] hig = statMax[kPop] end else: begin mid = stat50[kPop] low = stat25[kPop] hig = stat75[kPop] end endcase ;--- data cgPlot, /OverPlot, iArr, ssoHisto[kPop,*]*mid, color=popC[kPop], psym=10 ; cgErrPlot, iArr, ssoHisto[kPop,*]*low, ssoHisto[kPop,*]*hig, color=popC[kPop] cgPlot, /OverPlot, iArr, ssoHisto[kPop,*]*mid, color=popC[kPop], psym=popS[kPop] obsCal[kPop,*] = calHisto[kPop,*] * mid ;-- key xKey = 40. yKey = 1e5 xShi = 5 yKey = cgLogGen( nbPopPrint, start=1.25e4, finish=1.25e6) yKey2= cgLogGen( nbPopPrint, start=1e4, finish=1e6) xLen = 5 cgPlot, /OverPlot, xKey+[0,xLen], [yKey[kMain],yKey[kMain]], color=popC[kPop] cgPlot, /OverPlot, xKey+0.5*xLen, yKey[kMain], color=popC[kPop], psym=popS[kPop] cgText, xKey+xLen*1.2, yKey2[kMain], popName[kPop] if kMain eq 0 then begin cgErrPlot, xKey+0.5*xLen, 5e2*(hig/mid), 5e2*(mid/low) cgText, xKey+1.2*xLen, 2e3, 'Typical uncertainty' endif ; print, xKey+kMain*xShi, yKey*(hig/mid), yKey*(mid/low) ;hig/mid, mid/low print, popName[kPop], total(obsCal[kPop,*]), format='(A-15,1x,I7)' kMain++ endif endfor cgPlot, iArr, histoCal, /Normal, /NoErase, psym=10, color='Dark gray', $ xStyle=1, xRange=xRang, $ yStyle=1, yRange=[1, 1000], /yLog, $ yTitle='CalFields (#)', $ xTickInterval=5, xMinor=5, xTickName=replicate(' ',20), $ yTickInterval=100, yMinor=5, $ xTicklen=0.1, $ position=[0.09,0.80,0.96,0.98] cgPS_close, /PNG, delete_PS=0 cgPS_open, Filename=dirExten+'Euclid-Calib.eps', /metric, /decomposed, /encapsulated, $ xSize=35, ySize=24, language_level=2, /quiet xRang=[0,70] cgPlot, 0,0, /NoData, /Normal, $ xStyle=1, xRange=xRang, $ yStyle=1, yRange=[1, 1e6], /yLog, $ xTitle='Ecliptic latitude (!Uo!N)', $ yTitle='Cumulative number of SSOs in calibration fields', $ xTickInterval=5, xMinor=5, $ yTickInterval=1, yMinor=5, $ position=[0.09,0.105,0.96,0.96] for kPop=0, nbPop-1 do begin if popMain[kPop] eq 1 then begin cgPlot, /OverPlot, iArr, obsCal[kPop,*], color=popC[kPop], psym=10 cgPlot, /OverPlot, iArr, obsCal[kPop,*], color=popC[kPop], psym=popS[kPop] endif endfor cgPS_close, /PNG, delete_PS=0 ; forprint, iArr, (1.*histoI)*stat50[kPop]/totI, $ ; (1.*histoI)*stat25[kPop]/totI, (1.*histoI)*stat75[kPop]/totI, $ ; textout=dirExten+'euclid-latitude.csv', $ ; format='(I2,",",I7,",",I7,",",I7)', /Silent, $ ; comment='Latitude, #SSO, Q25, Q75' ; end