; ; statistics used for revision ; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;----What to do readDISCO = 01 readRS = 01 kPopStart = 0 debug = 0 ;- 0 Do not create the eps figure, but plot more stuff on screen ;- 1 check CSD for each pop ;- 2 check SFD before convolution with euclid ;- 3 check h2v & euclid detection PDF if debug ne 0 then set_plot,'X' ;--I.1-- Define Directories ------------------------------------------------------------------; spawn, 'hostname', host case strTrim(host) of 'hyperion': dirEuclid = '/data/euclid/' ; 'endymion': dirEuclid = '/home/bcarry/work/data/euclid/' 'endymion': dirEuclid = '/data/euclid/' else: stop endcase dirFilter = dirEuclid+'filters/' dirTaxo = dirEuclid+'taxonomy/' dirStat = dirEuclid + 'stat/' dirSFD = dirStat + 'SFD/' dirFrac = dirEuclid + 'fraction/' dirPop = dirEuclid + 'synthetic/' dirRS = dirEuclid + 'rs/' ;--I.2-- Broad Characteristics of the Survey -------------------------------------------------; sLimit = 17.0 ;-- Saturation vLimit = 24.5 ;-- YJH -- 24.5 in V, but smearing... binH = 0.25 ;- Mag binElong = 0.5 ;- Degree minH = 2 maxH = 30 ;?? hRang = [2,30] ;-SFD xRang = [2,24] ;-plot ;--I.3-- Solar Elongation of the Survey ------------------------------------------------------; if readRS eq 1 or not keyword_set(rs) then $ rs=euclidReadRS(dirRS+'wide.csv') cgHistoPlot, rs.SAA, binSize=binElong, loc=elong, histData=elongPDF euclidElong = [ [elong],[elongPDF] ] wdelete ;--I.4-- V-Mag Completeness Curves -----------------------------------------------------------; ;--I.4.1-- Mid-Range Completeness (Vis-Only in Kummel) readcol, dirStat+'completeness-Vis.csv',vMag,comp,delim=',', /Silent vInterp = 22 + binH*findgen( (27-22)/binH ) cInterp = interpol(comp, vMag, vInterp, /Spline) completeness = { V:vInterp, low:cInterp, mid:cInterp, high:cInterp } ;--I.4.2-- Low-Range Completeness (NISP-Only in Kummel) readcol, dirStat+'completeness-NISP.csv',vMag,comp,delim=',', /Silent vInterp = 22 + binH*findgen( (27-22)/binH ) cInterp = interpol(comp, vMag, vInterp, /Spline) completeness.low=cInterp ;--I.4.3-- High-Range Completeness (co-Added in Kummel) readcol, dirStat+'completeness-coAdd.csv',vMag,comp,delim=',', /Silent vInterp = 22 + binH*findgen( (27-22)/binH ) cInterp = interpol(comp, vMag, vInterp, /Spline) completeness.high=cInterp ; cgPlot, vMag,comp, xr=[22,27],yr=[0,1],psym='Filled Circle' ; cgPlot, /OverPlot, completeness.v, completeness.low, linestyle=1 ; cgPlot, /OverPlot, completeness.v, completeness.mid ; cgPlot, /OverPlot, completeness.v, completeness.high, linestyle=2 ;stop ; nea = ssoHtoV(pop='NEA',elong=euclidElong) ; mba = ssoHtoV(pop='MB',elong=euclidElong) ; mca = ssoHtoV(pop='MC',elong=euclidElong) ; ; cgPlot, mba.mag,mba.pdf ; cgPlot, nea.mag,nea.pdf,/OverPlot, color='gray' ; cgPlot, mca.mag,mca.pdf,/OverPlot, color='orange' if debug eq 3 then begin set_plot, 'x' cgplot, 0,0,xr=[0,27], yr=[0,1] compCurve = fltarr(nbH,3) for kC=0, nbComp-1 do begin hLimFaint = completeness.V[kC] - H2V.mag envFaint = interpol( total(H2V.pdf,/Cumulative), hLimFaint, sfdM.h ) compCurve[*,0] += envFaint * completeness.low[kC] compCurve[*,1] += envFaint * completeness.mid[kC] compCurve[*,2] += envFaint * completeness.high[kC] cgPlot, /OverPlot, sfdm.H, compCurve[*,1] endfor compCurve[where(compCurve gt 1)] =1. ;-Avoid probability > 1 ; compCurve[where(compCurve le 1e-3)]=0 ;-Clean long tails of probability cgPlot, /OverPlot, sfdm.H, compCurve[*,0],color='blue' cgPlot, /OverPlot, sfdm.H, compCurve[*,2],color='red' endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- 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) ; nbPop=1 ;--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 ;--III.4-- Output Array & Selected Populationss ----------------------------------------------; ecStat = fltarr(nbPop,2,3,3) ;-0 Population ;-1 To Be Observed/Discovered ;-2 Low/Mid/High Estimate of population ;-3 Low/Mid/High Estimate of completeness ecWhisk = fltarr(nbPop,5) popPrint = where( popMain eq 1, nbPopPrint ) print, 'Populations considered: ('+strtrim(string(nbPop),2)+'): ', popName ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- SSO CSD and Discovery -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Plot parameters --------------------------------------------------------------------; ;--IV.1.1-- Open PS File if debug eq 0 then $ cgPS_open, Filename=dirStat+'Euclid-Stat.eps', /metric, /decomposed, /encapsulated, $ xSize=35, ySize=40.5, language_level=2, /quiet ;--IV.1.2-- Graphic Regions yBot = 0.05 ySep = 0.03 yTop = 0.97 yLen = (yTop-yBot-ySep)/2.; + yBot yTop2= 0.977 ;--IV.1.3-- Location of Keyh Caption xKey = 3 xLen = 1 yKey = cgLogGen( 5, start=2.7e5, finish=2.7e7) yKey2= cgLogGen( 5, start=2e5, finish=2e7) kPrint=0 ;--IV.1.4-- Log -> Linear Y-position for Key yKeyLine = cgLogGen( 8, start=8.7e1, finish=4.7e4) yKeyLine2= cgLogGen( 8, start=8e1, finish=4e4) ;--IV.1.5-- Key for Discovery yBotDisco = 4 yTopDisco = 20 yTxtDisco = 12 xDisco = 19 ;--IV.1.6-- Create Lower Graphic cgPlot, 0,0, /NoData, /yLog, /Normal, $ xTitle='Absolute magnitude (H)', $ yTitle='Cumulative number N( 1 ; compCurve[where(compCurve le 1e-3)]=0 ;-Clean long tails of probability envFaint = compCurve if debug eq 3 then begin cgPlot, /OverPlot, sfdm.H, compCurve[*,0],color='blue' cgPlot, /OverPlot, sfdm.H, compCurve[*,2],color='red' stop endif ;--V.4-- Compute Detection Probability -----------------------------------------------------; ;--V.4.1-- Both Envelops envPDF = envFaint ;-Envelop PDF is the 3 faint version envPDF[*,0] *= envBright ;-Low faint completeness * bright end envPDF[*,1] *= envBright ;-Mid faint completeness * bright end envPDF[*,2] *= envBright ;-High faint completeness * bright end ;--V.4.2-- Whiskers for Bright End w0 = max( where( envBright le 0.01 )) trash = min( abs(envBright-0.25), w25 ) trash = min( abs(envBright-0.50), w50 ) trash = min( abs(envBright-0.75), w75 ) w100 = min( where(envBright eq 1 )) whiskerBright = [ w0,w25,w50,w75,w100 ] ;--V.4.3-- Whiskers for Faint End w0 = min( where( envFaint[*,1] le 0.01 )) trash = min( abs(envFaint[*,1]-0.25), w25 ) trash = min( abs(envFaint[*,1]-0.50), w50 ) trash = min( abs(envFaint[*,1]-0.75), w75 ) w100 = max( where(envFaint[*,1] eq 1 )) whiskerFaint = [w0,w25,w50,w75,w100 ] ecWhisk[kPop,*] = whiskerFaint ;--V.5-- Cumulative Distribution of Observation/Discoveries by Euclid ----------------------; ;--V.5.1-- Low Estimate cumulL_L = total( sfdL.histo.model * envPDF[*,0], /Cumulative) cumulL_M = total( sfdL.histo.model * envPDF[*,1], /Cumulative) cumulL_H = total( sfdL.histo.model * envPDF[*,2], /Cumulative) ecStat[kPop,0,0,0] = max(cumulL_L) ecStat[kPop,0,0,1] = max(cumulL_M) ecStat[kPop,0,0,2] = max(cumulL_H) ecStat[kPop,1,0,0] = max(cumulL_L) - max(sfdL.cumul.obs) ecStat[kPop,1,0,1] = max(cumulL_M) - max(sfdL.cumul.obs) ecStat[kPop,1,0,2] = max(cumulL_H) - max(sfdL.cumul.obs) ;--V.5.2-- Mid Estimate cumulM_L = total( sfdM.histo.model * envPDF[*,0], /Cumulative) cumulM_M = total( sfdM.histo.model * envPDF[*,1], /Cumulative) cumulM_H = total( sfdM.histo.model * envPDF[*,2], /Cumulative) ecStat[kPop,0,1,0] = max(cumulM_L) ecStat[kPop,0,1,1] = max(cumulM_M) ecStat[kPop,0,1,2] = max(cumulM_H) ecStat[kPop,1,1,0] = max(cumulM_L) - max(sfdM.cumul.obs) ecStat[kPop,1,1,1] = max(cumulM_M) - max(sfdM.cumul.obs) ecStat[kPop,1,1,2] = max(cumulM_H) - max(sfdM.cumul.obs) ;--V.5.3-- High Estimate cumulH_L = total( sfdH.histo.model * envPDF[*,0], /Cumulative) cumulH_M = total( sfdH.histo.model * envPDF[*,1], /Cumulative) cumulH_H = total( sfdH.histo.model * envPDF[*,2], /Cumulative) ecStat[kPop,0,2,0] = max(cumulH_L) ecStat[kPop,0,2,1] = max(cumulH_M) ecStat[kPop,0,2,2] = max(cumulH_H) ecStat[kPop,1,2,0] = max(cumulH_L) - max(sfdH.cumul.obs) ecStat[kPop,1,2,1] = max(cumulH_M) - max(sfdH.cumul.obs) ecStat[kPop,1,2,2] = max(cumulH_H) - max(sfdH.cumul.obs) print, popName[kPop] print, reform(ecStat[kPop,0,0,*]) print, reform(ecStat[kPop,0,1,*]) print, reform(ecStat[kPop,0,2,*]) print, '' print, reform(ecStat[kPop,1,0,*]) print, reform(ecStat[kPop,1,1,*]) print, reform(ecStat[kPop,1,2,*]) ;stop ; set_plot,'x' ; cgplot, sfdm.h,envPDF[*,0],xr=[0,27] ; cgplot, sfdm.h,envPDF[*,1], /over, color='gray' ; cgplot, sfdm.h,envPDF[*,2], /over, color='red' ;stop ;--V.6-- Plot Euclid Cumulative Distribution -----------------------------------------------; ; cgPlot, /OverPlot, sfdM.h, cumulM_L, color=popC[kPop], linestyle=3 ;-Very similar to Mid ; cgPlot, /OverPlot, sfdM.h, cumulM_H, color=popC[kPop], linestyle=3 ;-Very similar to Mid cgPlot, /OverPlot, sfdM.h, cumulM_M, color=popC[kPop], linestyle=3 cgPlot, /OverPlot, sfdM.h[whiskerFaint], cumulM_M[whiskerFaint], color=popC[kPop], psym='Open Circle' ; cgText, sfdM.h[whiskerFaint]+0.15, cumulM_M[whiskerFaint], color=popC[kPop], $ ; strTrim(string([0,25,50,75,100],format='(I)'),2), charSize=1 ; cgPlot, /OverPlot, sfdM.h[whiskerBright], cumulM_M[whiskerBright], color=popC[kPop], psym='Open Square' ; cgText, sfdM.h[whiskerBright]+0.15, cumulM_M[whiskerBright], color=popC[kPop], $ ; strTrim(string([0,25,50,75,100],format='(I)'),2), charSize=1 cgPlot, /OverPlot, sfdM.h[whiskerFaint[0]], sfdM.cumul.obs[whiskerFaint[0]], color=popC[kPop], psym='Filled Circle' cgPlot, /OverPlot, sfdM.h[whiskerFaint[[0,0]]], $ [sfdM.cumul.obs[whiskerFaint[0]],cumulM_M[whiskerFaint[0]]], color=popC[kPop] ;- ;--V.6-- Plot Euclid Cumulative Distribution -----------------------------------------------; ;- cgPlot, /OverPlot, sfdM.h, obsM.cumul.model, color=popC[kPop], linestyle=3 ;- ;- cgPlot, /OverPlot, sfdM.h[whiskerFaint], obsM.cumul.model[whiskerFaint], color=popC[kPop], psym='Open Circle' ;- ;- cgText, sfdM.h[whiskerFaint]+0.15, obsM.cumul.model[whiskerFaint], color=popC[kPop], $ ;- strTrim(string([0,25,50,75,100],format='(I)'),2), charSize=1 ;- ;- ;- cgPlot, /OverPlot, sfdM.h[whiskerBright], obsM.cumul.model[whiskerBright], color=popC[kPop], psym='Open Square' ;- ;- cgText, sfdM.h[whiskerBright]+0.15, obsM.cumul.model[whiskerBright], color=popC[kPop], $ ;- strTrim(string([0,25,50,75,100],format='(I)'),2), charSize=1 ;- ;- ;- ;- cgPlot, /OverPlot, sfdM.h[whiskerFaint[0]], obsM.cumul.obs[whiskerFaint[0]], color=popC[kPop], psym='Filled Circle' ;- ;- cgPlot, /OverPlot, sfdM.h[whiskerFaint[[0,0]]], $ ;- [obsM.cumul.obs[whiskerFaint[0]],obsM.cumul.model[whiskerFaint[0]]], color=popC[kPop] ;- ;- if popMain[kPop] eq 1 then kPrint++ ; if kpop eq 1 then goto, OUT endif ;-- end of loop over pop & test popmain OUT: 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 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- VI -- Disk I/O -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--VI.1-- Write Stats on Disk ----------------------------------------------------------------; forprint, popName, $ ecStat[*,0,0,0], ecStat[*,0,0,1], ecStat[*,0,0,2], $ ;-Obs / Lpop / L+M+H comp ecStat[*,0,1,0], ecStat[*,0,1,1], ecStat[*,0,1,2], $ ;-Obs / Mpop / L+M+H comp ecStat[*,0,2,0], ecStat[*,0,2,1], ecStat[*,0,2,2], $ ;-Obs / Hpop / L+M+H comp ecStat[*,1,0,0], ecStat[*,1,0,1], ecStat[*,1,0,2], $ ;-Disco / Lpop / L+M+H comp ecStat[*,1,1,0], ecStat[*,1,1,1], ecStat[*,1,1,2], $ ;-Disco / Mpop / L+M+H comp ecStat[*,1,2,0], ecStat[*,1,2,1], ecStat[*,1,2,2], $ ;-Disco / Hpop / L+M+H comp format='(A-20,18(", ",E10.3))', $ textout=dirStat+'euclid-stat.csv', /Silent,$ comment='Population, '+$ 'Obs_Lpop_L_comp, Obs_Lpop_M_comp, Obs_Lpop_H_comp, '+$ 'Obs_Mpop_L_comp, Obs_Mpop_M_comp, Obs_Mpop_H_comp, '+$ 'Obs_Hpop_L_comp, Obs_Hpop_M_comp, Obs_Hpop_H_comp, '+$ 'Disco_Lpop_L_comp, Disco_Lpop_M_comp, Disco_Lpop_H_comp, '+$ 'Disco_Mpop_L_comp, Disco_Mpop_M_comp, Disco_Mpop_H_comp, '+$ 'Disco_Hpop_L_comp, Disco_Hpop_M_comp, Disco_Hpop_H_comp, ' forprint, popName, $ sfdM.h[ecWhisk[*,0]], $ sfdM.h[ecWhisk[*,1]], $ sfdM.h[ecWhisk[*,2]], $ sfdM.h[ecWhisk[*,3]], $ sfdM.h[ecWhisk[*,4]], $ format='(A-20,5(", ",F5.2))', textout=dirStat+'euclid-faintH.csv', /Silent,$ comment='Population, H0, H25, H50, H75, H100' ; ;--VI.2-- Get Survey Fractions ---------------------------------------------------------------; ; readcol, dirFrac+'sso.fraction', delimiter=',', /Silent, $ ; fPop, fAcr, fWide, fCal, fDeep, fTot, format='(A,A,F,F,F,F)' ; stop ; 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, NallLow, NdiscoLow, NallMid, NdiscoMid, NallHigh, NdiscoHigh' ; 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 ;- code for getting quick & dirty mag of # of SSOs ; ; readcol,'sso.fraction',pop,aa,w,c,d,t,delim=',', format='(A,A,D,D,D,D)' ; forprint, pop, 100.*w/t,100.*c/t,1e6*d/t, 100*(w+d+c)/t,t/15,format='(A-14,2x,F5.1,2x,F5.1,2x,I2,2x,F5.1,2x,I6)' ; NEA 7.0 0.7 58 7.7 16061 ; MC 8.8 0.5 12 9.3 15487 ; MB 1.5 0.6 0 2.1 675046 ; Trojan 5.0 0.4 9 5.4 6761 ; Centaur 11.8 0.4 0 12.2 469 ; KBO 4.8 0.3 28 5.1 2330