;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-What to do checkFrac = 01 ;-Check fraction within Euclid RS ecExtent = 0 ;-What if Survey was extended toward the ecliptic? exOnlyFig = 0 ;-Only run extension figures ;---------- 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/' 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. ; rs=euclidReadRS( dirRS+'rs.csv' ) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- SSO populations -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--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) 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, oLow, dLow, oMid, dMid, oHig, dHig, format='(A,F,F,F,F,F,F)' statObs = oMid statDisco = dMid LowObs = oLow - oMid LowDisco = dLow - dMid HighObs = oHig - oMid HighDisco = dHig - dMid ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Euclid Reference Survey -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1 -- Survey definitions rsName = ['north', 'south'] rsNum = [1, 2] rs = replicate({id:'', $ ;-Region Id ns:0, $ ;-North=0 South=1 num:0, $ ;-North 1/2, South 1/2 lim:replicate({n:0, $ ;-Set to 1 if bin is in region ra:0., $ ;-RA of the bin minDec:0., $ ;-Lower limit in Dec maxDec:0.},$ ;-Higher limit in Dec 361)},4) ;--III.2 -- Loop over subset : north/south defined by hand in topcat kRS=0 for kName=0,1 do $ for kNum=0,1 do begin ;--III.2.1-- Read survey fileName = 'rs-'+rsName[kName]+string(rsNum[kNum],format='(I1)')+'.csv' rs[kRS].id = rsName[kName]+string(rsNum[kNum],format='(I1)') rs[kRS].ns = kName rs[kRS].num = kNum cur=euclidReadRS( dirRS+fileName ) ;--III.2.2-- Histogram in RA for dec limits (step = 1 deg) cgHistoPlot, cur.lon, binSize=1., loc=raArr, histData=nArr, minInput=0, maxInput=360 nRA = n_elements( raArr ) n1Place = where( nArr gt 0, nN ) rs[kRS].lim[n1Place].n = 1 ;--III.2.3-- Find min/max Dec lim = fltarr(nRA,2) for k=0, nRa-2 do begin range = where ( cur.lon ge raArr[k] and cur.lon le raArr[k+1] ) lim[k,0] = min( cur.lat[range], max=maxi ) lim[k,1] = maxi endfor ;--III.2.4-- Set min/max to poles if needed case kName of ;----- NORTH 0: begin if kNum eq 0 then lim[*,1] = 90. end ;----- SOUTH 1: begin if kNum eq 0 then lim[*,0] = -90. end endcase rs[kRS].lim.ra = raArr rs[kRS].lim.minDec = lim[*,0] rs[kRS].lim.maxDec = lim[*,1] ;--III.2.5-- Quality control cgplot, cur.lon, cur.lat, psym=4,xr=[0,360] cgPlot, rs[kRS].lim.ra, rs[kRS].lim.minDec, color='Red',/Over, psym=10, thick=3 cgPlot, rs[kRS].lim.ra, rs[kRS].lim.maxDec, color='Blue',/Over, psym=10, thick=3 ok = where( rs[kRS].lim.n eq 1, nOk ) if kName eq 0 and kNum eq 0 then pLat = 95 if kName eq 0 and kNum eq 1 then pLat = 10 if kName eq 1 and kNum eq 0 then pLat =-95 if kName eq 1 and kNum eq 1 then pLat =-15 cgplot, rs[kRS].lim[ok].ra, replicate(pLat,nOk ), psym=4, /over kRS++ endfor wdelete calib=euclidReadRS( dirRS+'calib.csv' ) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Determine reference epochs -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; readcol, dirEuclid+'time.sso', ep, format='(A)', /silent nbEp = n_elements( ep ) nbEp = 2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Check % of SSO in Euclid RS = f(pop,time) -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; IF checkFRAC eq 1 then begin ;--V.1.1-- Bookkeeping arrays stat = lonarr(nbPop, nbEp, 6) frac = fltarr(nbPop, nbEp, 6) ; tot = lonarr(nbPop, nbEp, 2) ; mainVec = intarr(nbPop) ssoNum = fltarr(nbPop,6,5) ;-Number of known object - in/out etc ssoFrac= fltarr(nbPop,6,5) ;-Fraction of known objects ssoTot = fltarr(nbPop,5) ;-Tot observed by euclid nbtot=1.d statCAL = lonarr(nbPop, nbEp, 6) fracCAL = fltarr(nbPop, nbEp, 6) calFrac= fltarr(nbPop,6,5) ;-Fraction of known objects IN THE CALIBRATION FIELDS ;--V.1.2-- Result arrays ; openw, idEx, dirFrac+'tab-frac.tex', /Get_Lun ;--V.2-- Loop over populations for kPop=0, nbPop-1 do begin ;--V.2.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' 'Detached': suffix = 'kbo-detached' 'SDO': suffix = 'kbo-sdo' 'Resonant': suffix = 'kbo-res' 'Inner Classical Belt': suffix = 'kbo-icb' 'Main Classical Belt': suffix = 'kbo-mcb' 'Outer Classical Belt': suffix = 'kbo-ocb' 'Comet': suffix = 'com' endcase ;--V.3-- Loop over epochs for kEp=0, nbEp-1 do begin ;--V.3.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 nbSSO = n_elements(ra1) ;--V.3.2-- Blunt cuts on galactic / ecliptic latitude nbSSO = n_elements(id1) noGal = where( abs(gLat1) ge limGal, nbGal ) noEcl = where( abs(eLat1) ge limEcl, nbEcl ) surv = where( abs(gLat1) ge limGal and abs(eLat1) ge limEcl, nbSurv, ncomp=nbCom ) stat[kPop,kEp,*] = [nbSSO, nbGal, nbEcl, nbSurv, 0, nbCom] frac[kPop,kEp,*] = [nbSSO, nbGal, nbEcl, nbSurv, 0, nbCom]*1./nbSSO statCAL[kPop,kEp,*] = [nbSSO, nbGal, nbEcl, nbSurv, 0, nbCom] fracCAL[kPop,kEp,*] = [nbSSO, nbGal, nbEcl, nbSurv, 0, nbCom]*1./nbSSO ;--V.4-- Real extent of Euclid RS -------------------------------------------------- nbQ = 0L for kRS=0, 3 do begin ; if kRS le 1 then yr=[0,90] else yr=[-90,0] ; cgplot, rs[kRS].lim.ra, rs[kRS].lim.minDec,xr=[0,360], yr=yr, yst=1 ; cgplot, /Over, rs[kRS].lim.ra, rs[kRS].lim.maxDec ; cgPlot, /OverPlot, eLong1, eLat1, psym=3, color='Navy' ;--V.4.1-- For each RS subset (North1-2 + South1-2) count=0L locQ = where( rs[kRS].lim.n eq 1, nbLoc ) for kLoc=0, nbLoc-1 do begin locSSo = where( eLong1 ge rs[kRS].lim[locQ[kLoc]].ra and eLong1 le rs[kRS].lim[locQ[kLoc]].ra+1 and $ eLat1 ge rs[kRS].lim[locQ[kLoc]].minDec and eLat1 le rs[kRS].lim[locQ[kLoc]].maxDec, nbIn ) count+=nbIn ; cgplot, /over, eLong1[locSSO], eLat1[locSSO], psym=4, color='red' endfor nbQ += count endfor stat[kPop,kEp,4] = nbQ frac[kPop,kEp,4] = nbQ*1./nbSSO ;--V.5-- CALIBRATION FIELDS ------------------------------------------------------- nbCal = n_elements( calib.lon ) nbQ = 0L for kC=0, nbCal-1 do begin gCirc, 2, calib.lon[kC], calib.lat[kC], elong1, elat1, angDist angDist *= 3600. inFoV = where( angDist le 0.367423, nbInFoV ) nbQ += nbInFoV endfor print, ' ---------------[ ',nbQ, ' ] ------------------' statCAL[kPop,kEp,4] = nbQ fracCAL[kPop,kEp,4] = nbQ*1./nbSSO endfor ;--V.5-- Average Epoch for kVal=0,5 do begin qStat = dataBox( stat[kPop,*,kVal] ) qFrac = dataBox( frac[kPop,*,kVal] ) med = mean( frac[kPop,*,kVal] ) std = stddev(frac[kPop,*,kVal] ) ssoNum[kPop,kVal,*] = qStat ssoFrac[kPop,kVal,*] = qFrac ssoFrac[kPop,kVal,2] = med ssoFrac[kPop,kVal,1] = med-std ssoFrac[kPop,kVal,3] = med+std qStat = dataBox( statCAL[kPop,*,kVal] ) qFrac = dataBox( fracCAL[kPop,*,kVal] ) med = mean( fracCAL[kPop,*,kVal] ) std = stddev(fracCAL[kPop,*,kVal] ) calFrac[kPop,kVal,*] = qFrac calFrac[kPop,kVal,2] = med calFrac[kPop,kVal,1] = med-std calFrac[kPop,kVal,3] = med+std endfor ;frac = ssoFrac[kPop,4,*] ;------- Compute the total amount of target to be observed by Euclid ssoTot[kPop,*] = 1. * statObs[kPop] * ssoFrac[kPop,4,*] nbTot += ssoTot[kPop,2] head='\hspace{1em}' if popName[kPop] eq 'NEA' or $ popName[kPop] eq 'MC' or $ popName[kPop] eq 'MB' or $ popName[kPop] eq 'Trojan' or $ popName[kPop] eq 'Centaur' or $ popName[kPop] eq 'KBO' or $ popName[kPop] eq 'Comet' then begin head='' endif fracV= 100.*ssoFrac[kPop,4,2] fracL= fracV - 100.*ssoFrac[kPop,4,1] fracH= 100.*ssoFrac[kPop,4,3] - fracV print, '' print, calFrac[kPop,4,1:3] print, '' calV= 100.*calFrac[kPop,4,2] calL= calV - 100.*calFrac[kPop,4,1] calH= 100.*calFrac[kPop,4,3] - calV print, head+popName[kPop], ssoNum[kPop,0,2], $ toSciNote( statObs[kPop], lowObs[kPop], highObs[kPop] , digit=1), $ toSciNote( statDisco[kPop], lowDisco[kPop], highDisco[kPop], digit=1), $ toSciNote( fracV, fracL, fracH, power=0, digit=1), $ toSciNote( calV, calL, calH, power=0, digit=1), $ toSciNote( statObs[kPop] *ssoFrac[kPop,4,2], lowObs[kPop] *ssoFrac[kPop,4,2], highObs[kPop] *ssoFrac[kPop,4,2], digit=1), $ toSciNote( statDisco[kPop]*ssoFrac[kPop,4,2], lowDisco[kPop]*ssoFrac[kPop,4,2], highDisco[kPop]*ssoFrac[kPop,4,2], digit=1), $ format='(A-25," & ",I6," & $",'+$ 'A-40," $&$ ",'+$ 'A-40," $&$ ",'+$ 'A-20," $&$ ",'+$ 'A-20," $&$ ",'+$ 'A-40," $&$ ",'+$ 'A-40," $\\ ")' endfor order = [7,6,5,4,2,1,0] openw, id, dirFrac+'euclid-fraction.tex', /get_Lun for kPop=0, 6 do begin kP = order[kPop] printf, id, head+popName[kPop], ssoNum[kPop,0,2], $ toSciNote( statObs[kPop], lowObs[kPop], highObs[kPop] , digit=1), $ toSciNote( statDisco[kPop], lowDisco[kPop], highDisco[kPop], digit=1), $ toSciNote( fracV, fracL, fracH, power=0, digit=1), $ toSciNote( statObs[kPop] *ssoFrac[kPop,4,2], lowObs[kPop] *ssoFrac[kPop,4,2], highObs[kPop] *ssoFrac[kPop,4,2], digit=1), $ toSciNote( statDisco[kPop]*ssoFrac[kPop,4,2], lowDisco[kPop]*ssoFrac[kPop,4,2], highDisco[kPop]*ssoFrac[kPop,4,2], digit=1), $ format='(A-25," & ",I6," & $",'+$ 'A-40," $&$ ",'+$ 'A-40," $&$ ",'+$ 'A-20," $&$ ",'+$ 'A-40," $&$ ",'+$ 'A-40," $\\ ")' endfor ff= weighted_mean( ssoFrac[[order],4,2], statObs[order]) fracV= weighted_mean( 100.*ssoFrac[order,4,2] , statObs[order] ) fracL= weighted_mean( fracV - 100.*ssoFrac[order,4,1], statObs[order] ) fracH= weighted_mean( 100.*ssoFrac[order,4,3] - fracV, statObs[order] ) printf,id, ' \hline' printf,id, 'Total', total(ssoNum[order,0,2]), $ toSciNote( total(statObs[order]), total(lowObs[order] ), total(highObs[order] ), digit=1), $ toSciNote( total(statDisco[order]), total(lowDisco[order]), total(highDisco[order]), digit=1), $ toSciNote( fracV, fracL, fracH, power=0, digit=1), $ toSciNote( total(statObs[order] *ssoFrac[order,4,2]), total(lowObs[order] *ssoFrac[order,4,2]), total(highObs[order] *ssoFrac[order,4,2]), digit=1), $ toSciNote( total(statDisco[order]*ssoFrac[order,4,2]), total(lowDisco[order]*ssoFrac[order,4,2]), total(highDisco[order]*ssoFrac[order,4,2]), digit=1), $ format='(A-25," & ",I6," & $",'+$ 'A-40," $&$ ",'+$ 'A-40," $&$ ",'+$ 'A-20," $&$ ",'+$ 'A-40," $&$ ",'+$ 'A-40," $\\ ")' close, id help, id endif ;----------------------------------------------------------------------------------------------------; ;----------------------------------------------------------------------------------------------------; ;----------------------------------------------------------------------------------------------------; ;----------------------------------------------------------------------------------------------------; ;----------------------------------------------------------------------------------------------------; ;----------------------------------------------------------------------------------------------------; ;----------------------------------------------------------------------------------------------------; ;----------------------------------------------------------------------------------------------------; ;----------------------------------------------------------------------------------------------------; IF ecExtent eq 1 then begin if exOnlyFig ne 1 then begin nbEc = round(maxEc)+1 iLim = maxEc - indgen(nbEc) stat = lonarr(nbPop+1,nbEc,5) frac = fltarr(nbPop+1,nbEc,5) for kPop=0, nbPop-1 do begin readcol, dirSSO+suffix[kPop]+'-'+ep[0]+'.csv', delimiter=',', /Silent, $ id1, ra1, dec1, d1, mag1, h1, g1, class1, name1, glon1, glat1, elong1, elat1, $ format='(I,F,F,F,F,F,F,A,A,F,F,F,F)' nbSSO = n_elements(id1) for kEc=0, nbEc-1 do begin noGal = where( abs(gLat1) ge limGal, nbGal ) noEcl = where( abs(eLat1) ge iLim[kEc], nbEcl ) surv = where( abs(gLat1) ge limGal and abs(eLat1) ge iLim[kEc], nbSurv, ncomp=nbCom ) stat[kPop,kEc,*] = [nbSSO, nbGal, nbEcl, nbSurv, nbcom] frac[kPop,kEc,*] = [nbSSO, nbGal, nbEcl, nbSurv, nbcom]*1./nbSSO stat[nbPop,kEc,*] += [nbSSO, nbGal, nbEcl, nbSurv, nbcom] endfor forprint, iLim, stat[kPop,*,3], 100.*frac[kPop,*,3], $ stat[kPop,*,2], 100.*frac[kPop,*,2], $ format='(I2,2x,I7,2x,F5.2,2x,I7,2x,F7.2)', textout=dirExten+popName[kPop]+'.dat', $ /Silent, /NoComment endfor for kEc=0, nbEc-1 do $ frac[nbPop,kEc,*] = 1.*stat[nbPop,kEc,*]/stat[nbPop,kEc,0] endif cgPS_open, Filename=dirExten+'Euclid-Exten.eps', /metric, /decomposed, /encapsulated, $ xSize=35, ySize=18.9, language_level=2, /quiet xRang=[0,25] cgPlot, 0,0, /NoData, /Normal, $ xStyle=1, xRange=xRang, $ yStyle=1, yRange=[0, 100], $ xTitle='Ecliptic latitude threshold (!Uo!N)', $ yTitle='Fraction of population (%)', $ xTickInterval=5, xMinor=5, $ yTickInterval=20, yMinor=4, $ position=[0.08,0.105,0.96,0.98] xKey = 18 xLen = 1.5 yKey = 65 + indgen(nbPop)*5 fracSel = 3 pEuclid=where( iLim eq limEcl) pMax = max(100*frac[*,pEuclid,fracSel]) ; cgPlot, /OverPlot, [limEcl,limEcl], [pMax+3, 80], color='Gray', linestyle=2 cgPlot, /OverPlot, [limEcl,limEcl], [3, 97], color='Gray', linestyle=2 cgText, limEcl-0.1, pMax+5, 'Limit of Euclid survey', color='Gray', $ orient=90 for kPop=0, nbPop-1 do begin cgPlot, /OverPlot, iLim, 100*frac[kPop,*,fracSel], color=popC[kPop] cgPlot, /OverPlot, iLim, 100*frac[kPop,*,fracSel], color=popC[kPop], psym=popS[kPop] cgPlot, /OverPlot, xKey+[0,xLen], yKey[kPop]+[0,0], color=popC[kPop] cgPlot, /OverPlot, xKey+0.5*xLen, yKey[kPop], color=popC[kPop], psym=popS[kPop] cgText, xKey+1.1*xLen, yKey[kPop]-2, popName[kPop] endfor ;---- WHOLE population of course follows MBA ; cgPlot, /OverPlot, iLim, 100*frac[nbPop,*,fracSel], color='Dark gray' xNum = [5, 10, 15] xNum = [0, 3, 6, 9, 12, 15] for kN=0, n_elements(xNum)-1 do begin pNum = where( iLim eq xNum[kN]) numObs = frac[nbPop,pNum,fracSel]*total(stat50);totalNumber numE = strtrim(string(numObs, format='(E7.0)'),2) numN = strmid(numE,0,1) expN = strmid(numE,5,1) cgText, xNum[kN]-0.1, 95, numN+'.10!U'+expN, align=1, orient=90 endfor cgPS_close, /PNG, delete_PS=0 endif ; end