; ; statistics used for revision ; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;----What to do readDISCO = 0 readRS = 0 ;---------- 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/' dirFrac = dirEuclid + 'fraction/' dirPop = dirEuclid + 'synthetic/' ; dirStat = dirEuclid + 'gaia-stat/' ;--- for Gaia only ;--- survey limits limEcl = 15 limGal = 20. sLimit = 17 ;-Saturation vLimit = 24.5 ;-- YJH -- 24.5 in V, but smearing... binH = 0.25 ;-geometry of obsevration elongation=90 elongMin=87 elongMax=110 now = 2016 jdNow = date_conv( string(now,format='(I4)')+'-01-01', 'JULIAN' ) nbY = 10 hFrac=[1.00,0.75,0.50,0.25,0.00] lFrac=reverse(hFrac) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- 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) ;-- pop, obs/disco, Low/Mid/High ecStat = fltarr(nbPop,2,3) popPrint = where( popMain eq 1, nbPopPrint ) print, 'Populations considered: ('+strtrim(string(nbPop),2)+'): ', popName ;--III.2-- Read SSO discovery archive if readDISCO eq 1 or not keyword_set(disco) then $ disco=ssoDiscovery_read() res = where( strCmp( strMid(disco.orbit.sub,0,8), 'resonant', /fold ), nbRes ) disco[res].orbit.sub = 'Resonant' ;--III.3-- Read LSST discoveries readcol, dirStat+'lsst', lsstPop, lsstTot, format='(A,F)', /Silent lsstPop=strtrim(lsstPop,2) lsstY1 = 2021 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- SSO CSD and discovery -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Plot parameters cgPS_open, Filename=dirStat+'Euclid-Stat.eps', /metric, /decomposed, /encapsulated, $ xSize=35, ySize=40.5, language_level=2, /quiet yBot = 0.05 ySep = 0.03 yTop = 0.97 yLen = (yTop-yBot-ySep)/2.; + yBot yTop2= 0.977 xKey = 3 xLen = 1 yKey = cgLogGen( 5, start=2.7e5, finish=2.7e7) yKey2= cgLogGen( 5, start=2e5, finish=2e7) kPrint=0 yKeyLine = cgLogGen( 8, start=8.7e1, finish=4.7e4) yKeyLine2= cgLogGen( 8, start=8e1, finish=4e4) yBotDisco = 4 yTopDisco = 20 yTxtDisco = 12 xDisco = 19 xRang = [2,24] cgPlot, 0,0, /NoData, /yLog, /Normal, $ xTitle='Absolute magnitude (H)', $ yTitle='Cumulative number N( 80ish ', comModelM.cumul.model[pp] sfdM = comModelM sfdL = comModelL sfdH = comModelH sfdM.cumul.obs=interp hLow=15.5 low = where( sfdM.h le hLow ) sfdM.cumul.model[low] = sfdM.cumul.obs[low] low = where( sfdL.h le hLow ) sfdL.cumul.model[low] = sfdM.cumul.obs[low] low = where( sfdH.h le hLow ) sfdH.cumul.model[low] = sfdM.cumul.obs[low] ; cgPlot, sfdM.h, comModelM.cumul.model, /over, color='gray' ; cgPlot, sfdL.h, sfdL.cumul.model, /over, color='gray', lineStyle=2 ; cgPlot, sfdM.h, sfdM.cumul.model, /over, color='gray' ; cgPlot, sfdH.h, sfdH.cumul.model, /over, color='gray', lineStyle=2 end ;--IV.2.2/X-- Default: from ASTORB and PowerLaw else: begin sel = where( strcmp(disco.orbit.main,popName[kPop]) or strcmp(disco.orbit.sub,popName[kPop]),nbSel ) ref=[popHmin[kPop],popHmax[kPop]] sfdM = ssoSynthSFD(disco[sel].h, reference=ref, binSize=binH, range=[0,24]) saveSFD = sfdM sfdL=sfdM sfdH=sfdM end endcase if popMain[kPop] eq 1 then begin ;--IV.2.4-- OverPlot synthetic distribution cgPlot, /OverPlot, sfdM.h, sfdM.cumul.model, color=popC[kPop], linestyle=2 cgPlot, /OverPlot, sfdM.h, sfdM.cumul.obs, color=popC[kPop] ;--IV.2.5-- Compute H-V minH2V = ssoHtoV(pop=popName[kpop], elong=elongMin) avgH2V = ssoHtoV(pop=popName[kpop], elong=elongation) maxH2V = ssoHtoV(pop=popName[kpop], elong=elongMax) joint= [minH2V, avgH2V, maxH2V] ;-IV.2.6-- Saturation and Limiting magnitudes hLimit = round( (vLimit-joint)/binH) * binH goodH = where( hLimit le 30 and finite(hLimit) and hLimit ge popHmax[kPop] ) hLimit = hLimit[goodH] lLimit = round( (sLimit-joint)/binH) * binH goodL = where( lLimit le 30 and finite(lLimit) and lLimit ge 0.1 ) lLimit = lLimit[goodL] hBox = dataBox( hLimit) lBox = dataBox( lLimit) ecSFD = sfdM.cumul.model * 0. ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Compute Average Numbers -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; curSFD = sfdM ;--V.1-- Rising from Saturation ----------------------------------------------------------; for kQ=0,4 do begin pL = where( curSFD.h eq lBox[kQ] ) pL = pL[0] xQ = curSFD.h[pL] if kQ eq 0 then begin yRef = curSFD.cumul.obs[pL] xB = xQ ecSFD[pL]=0 endif else begin subL = where( curSFD.h gt xB and curSFD.h le xQ, nbSub) iFrac = interpol( lFrac[[kQ-1,kQ]], [xB,xQ], curSFD.h[subL] ) for kSub=0, nbSub-1 do begin if popName[kPop] eq 'Centaur' and curSFD.h[subL[kSub]] ge 9.5 then begin ;and xB ge 9.5 then begin yNew = curSFD.cumul.model[subL[kSub]] endif else begin yNew = curSFD.cumul.obs[subL[kSub]] endelse dY = (yNew - yRef)*iFrac[kSub] yRef += dY ecSFD[subL[kSub]] = dY endfor xB = xQ endelse endfor ;--V.2-- Follow Model SFD ----------------------------------------------------------------; pH = where( curSFD.h gt lBox[4] and curSFD.h le hBox[0], nbH ) for kH=0, nbH-1 do begin yNew = curSFD.cumul.model[pH[kH]] dY = (yNew - yRef) yRef += dY ecSFD[pH[kH]] = dY endfor ;--V.3-- Dropping from Limiting Magnitude ------------------------------------------------; for kQ=0,4 do begin pH = where( curSFD.h eq hBox[kQ] ) pH = pH[0] xQ = curSFD.h[pH] if kQ eq 0 then begin yRef = curSFD.cumul.model[pH] xB = xQ endif else begin subH = where( curSFD.h gt xB and curSFD.h le xQ, nbSub) iFrac = interpol( hFrac[[kQ-1,kQ]], [xB,xQ], curSFD.h[subH] ) for kSub=0, nbSub-1 do begin yNew = curSFD.cumul.model[subH[kSub]] dY = (yNew - yRef)*iFrac[kSub] yRef += dY ecSFD[subH[kSub]] = dY endfor xB = xQ endelse endfor ;--V.4-- Cumulative SFD Observed by Euclid -----------------------------------------------; ecCum = total(ecSFD,/cum) ecStat[kPop,0,1] = max(ecCum) ecStat[kPop,1,1] = max(ecCum)-max(curSFD.cumul.obs) ;--V.5-- Plotting - Middle SFD only ------------------------------------------------------; ;--V.5.1-- Clean strange points case popName[kPop] of 'MB': begin pH=where(curSFD.h eq 10.7500) ecCum[pH]=1 end else: endcase cgPlot, curSFD.h, ecCum, /OverPlot, lineStyle=3, color=popC[kPop] ;--V.5.2-- Limiting Quartiles and margin for Discovery for kQ=0,4 do begin pH = where( curSFD.h eq hBox[kQ] ) pH = pH[0] xQ = curSFD.h[pH] yQ = ecCum[pH] cgPlot, /OverPlot, xQ, yQ, psym='Open Circle', color=popC[kPop] if kQ eq 4 then begin yN = curSFD.cumul.obs[pH] cgPlot, /OverPlot, xQ, yN, psym='Filled Circle', color=popC[kPop] cgPlot, /OverPlot, [xQ,xQ], [yN,yQ], color=popC[kPop] endif endfor ;--V.5.2-- Limiting Quartiles and margin for Saturation for kQ=0,4 do begin pL = where( curSFD.h eq lBox[kQ] ) pL = pL[0] xQ = curSFD.h[pL] yQ = ecCum[pL] cgPlot, /OverPlot, xQ, yQ, psym='Open square', color=popC[kPop] endfor ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- VI -- Compute Lower Numbers -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; curSFD = sfdL ;--VI.1-- Rising from Saturation ---------------------------------------------------------; for kQ=0,4 do begin pL = where( curSFD.h eq lBox[kQ] ) pL = pL[0] xQ = curSFD.h[pL] if kQ eq 0 then begin yRef = curSFD.cumul.obs[pL] xB = xQ ecSFD[pL]=0 endif else begin subL = where( curSFD.h gt xB and curSFD.h le xQ, nbSub) iFrac = interpol( lFrac[[kQ-1,kQ]], [xB,xQ], curSFD.h[subL] ) for kSub=0, nbSub-1 do begin if popName[kPop] eq 'Centaur' and curSFD.h[subL[kSub]] ge 9.5 then begin ;and xB ge 9.5 then begin yNew = curSFD.cumul.model[subL[kSub]] endif else begin yNew = curSFD.cumul.obs[subL[kSub]] endelse dY = (yNew - yRef)*iFrac[kSub] yRef += dY ecSFD[subL[kSub]] = dY endfor xB = xQ endelse endfor ;--VI.2-- Follow Model SFD ---------------------------------------------------------------; pH = where( curSFD.h gt lBox[4] and curSFD.h le hBox[0], nbH ) for kH=0, nbH-1 do begin yNew = curSFD.cumul.model[pH[kH]] dY = (yNew - yRef) yRef += dY ecSFD[pH[kH]] = dY endfor ;--VI.3-- Dropping from Limiting Magnitude -----------------------------------------------; for kQ=0,4 do begin pH = where( curSFD.h eq hBox[kQ] ) pH = pH[0] xQ = curSFD.h[pH] if kQ eq 0 then begin yRef = curSFD.cumul.model[pH] xB = xQ endif else begin subH = where( curSFD.h gt xB and curSFD.h le xQ, nbSub) iFrac = interpol( hFrac[[kQ-1,kQ]], [xB,xQ], curSFD.h[subH] ) for kSub=0, nbSub-1 do begin yNew = curSFD.cumul.model[subH[kSub]] dY = (yNew - yRef)*iFrac[kSub] yRef += dY ecSFD[subH[kSub]] = dY endfor xB = xQ endelse endfor ;--VI.4-- Cumulative SFD Observed by Euclid ----------------------------------------------; ecCum = total(ecSFD,/cum) ecStat[kPop,0,0] = max(ecCum) ecStat[kPop,1,0] = max(ecCum)-max(curSFD.cumul.obs) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- VII -- Compute Higher Numbers -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; curSFD = sfdH ;--VII.1-- Rising from Saturation --------------------------------------------------------; for kQ=0,4 do begin pL = where( curSFD.h eq lBox[kQ] ) pL = pL[0] xQ = curSFD.h[pL] if kQ eq 0 then begin yRef = curSFD.cumul.obs[pL] xB = xQ ecSFD[pL]=0 endif else begin subL = where( curSFD.h gt xB and curSFD.h le xQ, nbSub) iFrac = interpol( lFrac[[kQ-1,kQ]], [xB,xQ], curSFD.h[subL] ) for kSub=0, nbSub-1 do begin if popName[kPop] eq 'Centaur' and curSFD.h[subL[kSub]] ge 9.5 then begin ;and xB ge 9.5 then begin yNew = curSFD.cumul.model[subL[kSub]] endif else begin yNew = curSFD.cumul.obs[subL[kSub]] endelse dY = (yNew - yRef)*iFrac[kSub] yRef += dY ecSFD[subL[kSub]] = dY endfor xB = xQ endelse endfor ;--VII.2-- Follow Model SFD --------------------------------------------------------------; pH = where( curSFD.h gt lBox[4] and curSFD.h le hBox[0], nbH ) for kH=0, nbH-1 do begin yNew = curSFD.cumul.model[pH[kH]] dY = (yNew - yRef) yRef += dY ecSFD[pH[kH]] = dY endfor ;--VII.3-- Dropping from Limiting Magnitude ----------------------------------------------; for kQ=0,4 do begin pH = where( curSFD.h eq hBox[kQ] ) pH = pH[0] xQ = curSFD.h[pH] if kQ eq 0 then begin yRef = curSFD.cumul.model[pH] xB = xQ endif else begin subH = where( curSFD.h gt xB and curSFD.h le xQ, nbSub) iFrac = interpol( hFrac[[kQ-1,kQ]], [xB,xQ], curSFD.h[subH] ) for kSub=0, nbSub-1 do begin yNew = curSFD.cumul.model[subH[kSub]] dY = (yNew - yRef)*iFrac[kSub] yRef += dY ecSFD[subH[kSub]] = dY endfor xB = xQ endelse endfor ;--VII.4-- Cumulative SFD Observed by Euclid ---------------------------------------------; ecCum = total(ecSFD,/cum) ecStat[kPop,0,2] = max(ecCum) ecStat[kPop,1,2] = max(ecCum)-max(curSFD.cumul.obs) endif if popMain[kPop] eq 1 then kPrint++ endfor cgText, xKey+xLen*0.5, yKey2[kPrint], 'Inner Solar System' ;--- Diameter range dRang = (1329/sqrt(0.20)) * 10^(-0.2 * xRang) cgAxis, xAxis=1, xRange=dRang, /xLog, xStyle=1, $ xTickName=[' ','100','10','1','0.1'] cgText, /Normal, 0.95, yTop2, 'km' cgText, /Normal, 0.05, yTop2, 'Diameter' kPrint=7 cgText, xKey+xLen*0.5, yKeyLine2[kPrint], 'Populations', color='Dark Gray' kPrint-- cgPlot, /OverPlot, xKey+[0,xLen], yKeyLine[[kPrint,kPrint]], thick=6, color='Dark Gray' cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'Known', color='Dark Gray' kPrint-- cgPlot, /OverPlot, xKey+[0,xLen], yKeyLine[[kPrint,kPrint]], thick=6, lineStyle=2, color='Dark Gray' cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'Synthetic', color='Dark Gray' kPrint-- cgPlot, /OverPlot, xKey+[0,xLen], yKeyLine[[kPrint,kPrint]], thick=6, lineStyle=3, color='Dark Gray' cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'by Euclid', color='Dark Gray' kPrint-- kPrint-- cgText, xKey+xLen*0.5, yKeyLine2[kPrint], 'Detections', color='Dark Gray' kPrint-- cgPlot, /OverPlot, xKey+xLen/2., yKeyLine[[kPrint]], color='Dark Gray', psym='Open Circle' cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'Photon-limited', color='Dark Gray' kPrint-- cgPlot, /OverPlot, xKey+xLen/2., yKeyLine[[kPrint]], color='Dark Gray', psym='Open square' cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'Saturated', color='Dark Gray' cgText, xDisco+xLen*1.2, yTxtDisco, 'Margin for' cgText, xDisco+xLen*1.2, yBotDisco, 'discovery' cgPS_close, /PNG, delete_PS=0 ; forprint, popName, ecStat[*,0,1], ecStat[*,1,1], ecStat[*,0,1], ecStat[*,1,1], $ ; format='(A-20,", ",I8,", ",I8,", ",E9.3,", ",E9.3)', $ ; textout=dirStat+'euclid-stat.csv', /Silent, comment='Population, Nall, Ndisco, NallE, NdiscoE' forprint, popName, $ ecStat[*,0,0], ecStat[*,1,0], $ ;-low ecStat[*,0,1], ecStat[*,1,1], $ ;-mid ecStat[*,0,2], ecStat[*,1,2], $ ;-high format='(A-20,6(", ",E10.3))', $ textout=dirStat+'euclid-stat.csv', /Silent, comment='Population, Nall, Ndisco, NallE, NdiscoE' ; forprint, popName, fracI, $ ; discoStat[*,0,0], $ ; discoStat[*,0,1], $ ; discoStat[*,0,2], $ ; discoStat[*,0,3], $ ; discoStat[*,0,4], $ ; format='(A-20,",",F6.4,",",E8.1,",",E8.1,",",E8.1,",",E8.1,",",E8.1)', $ ; textout=dirStat+'euclid-disco.csv', /Silent, $ ; comment='Population, f_i, Min, Q25, Q50, Q75, Max' ; end