;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-What to do checkFrac = 0 ;-Check fraction within Euclid RS ecExtent = 01 ;-What if Survey was extended toward the ecliptic? exOnlyFig = 0 ;-Only run extension figures ;---------- 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. ;Time = 2020-01-31T00:00:00 ;RA/DEC from Euclid at 90 Elong = 14 37 36.54012 -09 14 23.0000 ra90 = ten(14,37,36.5)*15. mod 180 ; 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, statObs, statDisco, format='(A,L,L)' ; statName, statF, statMin, stat25, stat50, stat75, statMax, format='(A,F,F,F,F,F,F)' ; totalNumber = total(nbAvail) ;--II.2-- 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)' ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- 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.ra, 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.ra ge raArr[k] and cur.ra le raArr[k+1] ) lim[k,0] = min( cur.dec[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.ra, cur.dec, 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 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- 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) ;--V.1.2-- Result arrays ;-old ssoNum = fltarr(nbPop,6,5) ;-Number of known object - in/out etc ;-old ssoFrac= fltarr(nbPop,6,5) ;-Fraction of known objects ;-old ssoTot = fltarr(nbPop,5) ;-Tot observed by euclid ssoNum = fltarr(nbPop,6) ;-Number of known object - in/out etc ssoFrac= fltarr(nbPop,6) ;-Fraction of known objects nbTot = 0. tNow = lonarr(nbPop) tSky = fltarr(nbPop,2,3) tEuc = fltarr(nbPop,2,3) tFS = fltarr(nbPop,3) 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 ;--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 endfor ;--V.5-- Average Epoch for kVal=0,5 do begin qStat = dataBox( stat[kPop,*,kVal] ) qFrac = dataBox( frac[kPop,*,kVal] ) ssoNum[kPop,kVal,*] = qStat ssoFrac[kPop,kVal,*] = qFrac endfor ;------- Compute the total amount of target to be observed by Euclid ;-simple median * median ssoTot[kPop,*] = 1. * statObs[kPop] * ssoFrac[kPop,4,*] ;-more details [25-50-75]fs # [25-50-75]avail vecAvail=[stat25[kPop], stat50[kPop], stat75[kPop]] vecFrac =ssoFrac[kPop,4,1:3] spread = vecAvail # vecFrac ssoTot[kPop,*] = dataBox( spread ) numE = strtrim(string(ssoTot[kPop,2], format='(E9.1)'),2) numN = float(strmid(numE,0,3)) expN = float(strmid(numE,5,3)) empty={val:'', num:0., exp:0., fac:0.} obsEC={q25:empty, q50:empty, q75:empty} str50 = strtrim(string(ssoTot[kPop,2], format='(E+9.2)'),2) num50 = float( strmid(str50,0,5) ) exp50 = float( strmid(str50,6,3) ) obsEC.q50.val = str50 obsEC.q50.num = num50 obsEC.q50.exp = exp50 str75 = strtrim(string(ssoTot[kPop,3]-ssoTot[kPop,2], format='(E+9.2)'),2) num75 = float( strmid(str75,0,5) ) exp75 = float( strmid(str75,6,3) ) fac75 = float(exp75)-float(exp50) obsEC.q75.val = str75 obsEC.q75.num = num75 obsEC.q75.exp = exp75 obsEC.q75.fac = fac75 str25 = strtrim(string(ssoTot[kPop,1]-ssoTot[kPop,2], format='(E+9.2)'),2) num25 = float( strmid(str25,0,5) ) exp25 = float( strmid(str25,6,3) ) fac25 = float(exp25)-float(exp50) obsEC.q25.val = str25 obsEC.q25.num = num25 obsEC.q25.exp = exp25 obsEC.q25.fac = fac25 tEuc[kPop,0,*]= [ssoTot[kPop,1], ssoTot[kPop,2], ssoTot[kPop,3]] ;------- Compute the total amount of DISCO by Euclid vecAvail=[disco25[kPop], disco50[kPop], disco75[kPop]] vecFrac =ssoFrac[kPop,4,1:3] spread = vecAvail # vecFrac ssoTot[kPop,*] = dataBox( spread ) numE = strtrim(string(ssoTot[kPop,2], format='(E9.1)'),2) numN = float(strmid(numE,0,3)) expN = float(strmid(numE,5,3)) empty={val:'', num:0., exp:0., fac:0.} discoEC={q25:empty, q50:empty, q75:empty} str50 = strtrim(string(ssoTot[kPop,2], format='(E+9.2)'),2) num50 = float( strmid(str50,0,5) ) exp50 = float( strmid(str50,6,3) ) discoEC.q50.val = str50 discoEC.q50.num = num50 discoEC.q50.exp = exp50 str75 = strtrim(string(ssoTot[kPop,3]-ssoTot[kPop,2], format='(E+9.2)'),2) num75 = float( strmid(str75,0,5) ) exp75 = float( strmid(str75,6,3) ) fac75 = float(exp75)-float(exp50) discoEC.q75.val = str75 discoEC.q75.num = num75 discoEC.q75.exp = exp75 discoEC.q75.fac = fac75 str25 = strtrim(string(ssoTot[kPop,1]-ssoTot[kPop,2], format='(E+9.2)'),2) num25 = float( strmid(str25,0,5) ) exp25 = float( strmid(str25,6,3) ) fac25 = float(exp25)-float(exp50) discoEC.q25.val = str25 discoEC.q25.num = num25 discoEC.q25.exp = exp25 discoEC.q25.fac = fac25 if discoEC.q25.num lt 0 then discoEC.q25.num=0 tEuc[kPop,1,*]= [ssoTot[kPop,1], ssoTot[kPop,2], ssoTot[kPop,3]] ;-------------- The total number for the whole sky empty={val:'', num:0., exp:0., fac:0.} obsSky={q25:empty, q50:empty, q75:empty} str50 = strtrim(string(stat50[kPop], format='(E+9.2)'),2) num50 = float( strmid(str50,0,5) ) exp50 = float( strmid(str50,6,3) ) obsSky.q50.val = str50 obsSky.q50.num = num50 obsSky.q50.exp = exp50 str75 = strtrim(string(stat75[kPop]-stat50[kPop], format='(E+9.2)'),2) num75 = float( strmid(str75,0,5) ) exp75 = float( strmid(str75,6,3) ) fac75 = float(exp75)-float(exp50) obsSky.q75.val = str75 obsSky.q75.num = num75 obsSky.q75.exp = exp75 obsSky.q75.fac = fac75 str25 = strtrim(string(stat25[kPop]-stat50[kPop], format='(E+9.2)'),2) num25 = float( strmid(str25,0,5) ) exp25 = float( strmid(str25,6,3) ) fac25 = float(exp25)-float(exp50) obsSky.q25.val = str25 obsSky.q25.num = num25 obsSky.q25.exp = exp25 obsSky.q25.fac = fac25 ;-------------- The total number for the whole sky TO DISCOVER!!! empty={val:'', num:0., exp:0., fac:0.} discoSky={q25:empty, q50:empty, q75:empty} str50 = strtrim(string(disco50[kPop], format='(E+9.2)'),2) num50 = float( strmid(str50,0,5) ) exp50 = float( strmid(str50,6,3) ) discoSky.q50.val = str50 discoSky.q50.num = num50 discoSky.q50.exp = exp50 str75 = strtrim(string(disco75[kPop]-disco50[kPop], format='(E+9.2)'),2) num75 = float( strmid(str75,0,5) ) exp75 = float( strmid(str75,6,3) ) fac75 = float(exp75)-float(exp50) discoSky.q75.val = str75 discoSky.q75.num = num75 discoSky.q75.exp = exp75 discoSky.q75.fac = fac75 str25 = strtrim(string(disco25[kPop]-disco50[kPop], format='(E+9.2)'),2) num25 = float( strmid(str25,0,5) ) exp25 = float( strmid(str25,6,3) ) fac25 = float(exp25)-float(exp50) discoSky.q25.val = str25 discoSky.q25.num = num25 discoSky.q25.exp = exp25 discoSky.q25.fac = fac25 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='' mainVec[kPop] = 1 endif print, head+popName[kPop], ssoNum[kPop,0,2], $ obsSky.q50.num, obsSky.q25.num * 10^obsSky.q25.fac, obsSky.q75.num * 10^obsSky.q75.fac, obsSky.q50.exp, $ discoSky.q50.num, discoSky.q25.num * 10^discoSky.q25.fac, discoSky.q75.num * 10^discoSky.q75.fac, discoSky.q50.exp, $ 100.*ssoFrac[kPop,4,2], $ 100.*( ssoFrac[kPop,4,1]-ssoFrac[kPop,4,2] ), $ 100.*( ssoFrac[kPop,4,3]-ssoFrac[kPop,4,2] ), $ obsEC.q50.num, obsEC.q25.num * 10^obsEC.q25.fac, obsEC.q75.num * 10^obsEC.q75.fac, obsEC.q50.exp, $ discoEC.q50.num, discoEC.q25.num * 10^discoEC.q25.fac, discoEC.q75.num * 10^discoEC.q75.fac, discoEC.q50.exp, $ format='(A-25," & ",I6," & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F6.1,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ \\")' printf, idEx, head+popName[kPop], ssoNum[kPop,0,2], $ obsSky.q50.num, obsSky.q25.num * 10^obsSky.q25.fac, obsSky.q75.num * 10^obsSky.q75.fac, obsSky.q50.exp, $ discoSky.q50.num, discoSky.q25.num * 10^discoSky.q25.fac, discoSky.q75.num * 10^discoSky.q75.fac, discoSky.q50.exp, $ 100.*ssoFrac[kPop,4,2], $ 100.*( ssoFrac[kPop,4,1]-ssoFrac[kPop,4,2] ), $ 100.*( ssoFrac[kPop,4,3]-ssoFrac[kPop,4,2] ), $ obsEC.q50.num, obsEC.q25.num * 10^obsEC.q25.fac, obsEC.q75.num * 10^obsEC.q75.fac, obsEC.q50.exp, $ discoEC.q50.num, discoEC.q25.num * 10^discoEC.q25.fac, discoEC.q75.num * 10^discoEC.q75.fac, discoEC.q50.exp, $ format='(A-27," & ",I6," & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F6.1,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ \\")' ;--- some bookeeping for total line tNow[kPop] = ssoNum[kPop,0,2] tSky[kPop,0,*]= [stat25[kPop],stat50[kPop],stat75[kPop]] tSky[kPop,1,*]= [disco25[kPop],disco50[kPop],disco75[kPop]] ; tEuc[kPop,0,*]= [ssoTot[kPop,1], ssoTot[kPop,2], ssoTot[kPop,3]] tFS [kPop,*] =[ 100.*ssoFrac[kPop,4,2], 100.*( ssoFrac[kPop,4,1]-ssoFrac[kPop,4,2] ), 100.*( ssoFrac[kPop,4,3]-ssoFrac[kPop,4,2] ) ] endfor print, ' \hline' index = where( mainVec eq 1, nbMain ) totNow = total( tNow[index] ) totSky = total( tSky[index,*,*], 1 ) totEuc = total( tEuc[index,*,*], 1 ) totFS = total( tFS[index,*], 1 ) / float(nbMain) ;-OBS SKY TOTAL str50 = strtrim(string(totSky[0,1], format='(E+9.2)'),2) num50 = float( strmid(str50,0,5) ) exp50 = float( strmid(str50,6,3) ) obsSky.q50.val = str50 obsSky.q50.num = num50 obsSky.q50.exp = exp50 str75 = strtrim(string(totSky[0,2]-totSky[0,1], format='(E+9.2)'),2) num75 = float( strmid(str75,0,5) ) exp75 = float( strmid(str75,6,3) ) fac75 = float(exp75)-float(exp50) obsSky.q75.val = str75 obsSky.q75.num = num75 obsSky.q75.exp = exp75 obsSky.q75.fac = fac75 str25 = strtrim(string(totSky[0,0]-totSky[0,1], format='(E+9.2)'),2) num25 = float( strmid(str25,0,5) ) exp25 = float( strmid(str25,6,3) ) fac25 = float(exp25)-float(exp50) obsSky.q25.val = str25 obsSky.q25.num = num25 obsSky.q25.exp = exp25 obsSky.q25.fac = fac25 ;-- DISCO SKY TOTAL str50 = strtrim(string(totSky[1,1], format='(E+9.2)'),2) num50 = float( strmid(str50,0,5) ) exp50 = float( strmid(str50,6,3) ) discoSky.q50.val = str50 discoSky.q50.num = num50 discoSky.q50.exp = exp50 str75 = strtrim(string(totSky[1,2]-totSky[1,1], format='(E+9.2)'),2) num75 = float( strmid(str75,0,5) ) exp75 = float( strmid(str75,6,3) ) fac75 = float(exp75)-float(exp50) discoSky.q75.val = str75 discoSky.q75.num = num75 discoSky.q75.exp = exp75 discoSky.q75.fac = fac75 str25 = strtrim(string(totSky[1,0]-totSky[1,1], format='(E+9.2)'),2) num25 = float( strmid(str25,0,5) ) exp25 = float( strmid(str25,6,3) ) fac25 = float(exp25)-float(exp50) discoSky.q25.val = str25 discoSky.q25.num = num25 discoSky.q25.exp = exp25 discoSky.q25.fac = fac25 ;-- OBS EUCLID TOTAL str50 = strtrim(string(totEuc[0,1], format='(E+9.2)'),2) num50 = float( strmid(str50,0,5) ) exp50 = float( strmid(str50,6,3) ) obsEC.q50.val = str50 obsEC.q50.num = num50 obsEC.q50.exp = exp50 str75 = strtrim(string(totEuc[0,2]-totEuc[0,1], format='(E+9.2)'),2) num75 = float( strmid(str75,0,5) ) exp75 = float( strmid(str75,6,3) ) fac75 = float(exp75)-float(exp50) obsEC.q75.val = str75 obsEC.q75.num = num75 obsEC.q75.exp = exp75 obsEC.q75.fac = fac75 str25 = strtrim(string(totEuc[0,0]-totEuc[0,1], format='(E+9.2)'),2) num25 = float( strmid(str25,0,5) ) exp25 = float( strmid(str25,6,3) ) fac25 = float(exp25)-float(exp50) obsEC.q25.val = str25 obsEC.q25.num = num25 obsEC.q25.exp = exp25 obsEC.q25.fac = fac25 ;-- DISCO EUCLID TOTAL str50 = strtrim(string(totEuc[1,1], format='(E+9.2)'),2) num50 = float( strmid(str50,0,5) ) exp50 = float( strmid(str50,6,3) ) discoEC.q50.val = str50 discoEC.q50.num = num50 discoEC.q50.exp = exp50 str75 = strtrim(string(totEuc[1,2]-totEuc[1,1], format='(E+9.2)'),2) num75 = float( strmid(str75,0,5) ) exp75 = float( strmid(str75,6,3) ) fac75 = float(exp75)-float(exp50) discoEC.q75.val = str75 discoEC.q75.num = num75 discoEC.q75.exp = exp75 discoEC.q75.fac = fac75 str25 = strtrim(string(totEuc[1,0]-totEuc[1,1], format='(E+9.2)'),2) num25 = float( strmid(str25,0,5) ) exp25 = float( strmid(str25,6,3) ) fac25 = float(exp25)-float(exp50) discoEC.q25.val = str25 discoEC.q25.num = num25 discoEC.q25.exp = exp25 discoEC.q25.fac = fac25 if discoEC.q25.num lt 0 then discoEC.q25.num=0 print, 'Total', totNow, $ obsSky.q50.num, obsSky.q25.num * 10^obsSky.q25.fac, obsSky.q75.num * 10^obsSky.q75.fac, obsSky.q50.exp, $ discoSky.q50.num, discoSky.q25.num * 10^discoSky.q25.fac, discoSky.q75.num * 10^discoSky.q75.fac, discoSky.q50.exp, $ 100.*totFS[1], $ 100.*( totFS[0]-totFS[1] ), $ 100.*( totFS[1]-totFS[1] ), $ obsEC.q50.num, obsEC.q25.num * 10^obsEC.q25.fac, obsEC.q75.num * 10^obsEC.q75.fac, obsEC.q50.exp, $ discoEC.q50.num, discoEC.q25.num * 10^discoEC.q25.fac, discoEC.q75.num * 10^discoEC.q75.fac, discoEC.q50.exp, $ format='(A-25," & ",I6," & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F6.1,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ \\")' printf, idEx, '\hline' printf, idEx, 'Total', totNow, $ obsSky.q50.num, obsSky.q25.num * 10^obsSky.q25.fac, obsSky.q75.num * 10^obsSky.q75.fac, obsSky.q50.exp, $ discoSky.q50.num, discoSky.q25.num * 10^discoSky.q25.fac, discoSky.q75.num * 10^discoSky.q75.fac, discoSky.q50.exp, $ 100.*totFS[1], $ 100.*( totFS[0]-totFS[1] ), $ 100.*( totFS[1]-totFS[1] ), $ obsEC.q50.num, obsEC.q25.num * 10^obsEC.q25.fac, obsEC.q75.num * 10^obsEC.q75.fac, obsEC.q50.exp, $ discoEC.q50.num, discoEC.q25.num * 10^discoEC.q25.fac, discoEC.q75.num * 10^discoEC.q75.fac, discoEC.q50.exp, $ format='(A-27," & ",I6," & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F6.1,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ \\")' free_lun, idEx forprint, popName, tNow, $ ;-Pop + #now tSky[*,0,0], tSky[*,0,1], tSky[*,0,2], $ ;-Q25-50-75: obs/sky tSky[*,1,0], tSky[*,1,1], tSky[*,1,2], $ ;-Q25-50-75: disco/sky tFS[*,0], tFS[*,1]+tFS[*,0], tFS[*,2]+tFS[*,0], $;-Q25-50-75: fraction tEuc[*,0,0], tEuc[*,0,1], tEuc[*,0,2], $ ;-Q25-50-75: obs/euclid tEuc[*,1,0], tEuc[*,1,1], tEuc[*,1,2], $ ;-Q25-50-75: disco/euclid textout=dirFrac+'euclid-fraction.csv', /Silent, $ format='(A-15,", ",I7,6(", ",E13.6),3(", ",F7.4),6(", ",E13.6))' ;---old numE = strtrim(string(ssoTot[0,2], format='(E9.1)'),2) ;---old numN = strmid(numE,0,3) ;---old expN = strmid(numE,5,3) ;---old ;---old totalStat = total( stat, 1 ) ;---old numForAll = dataBox( totalStat[*,0] ) ;---old insForAll = dataBox( totalStat[*,4] ) ;---old fracForAll = 100.*insForAll/numForAll ;---old ;---old ;---old ;---old tot25 =total(stat25[index]) ;---old tot50 =total(stat50[index]) ;---old tot75 =total(stat75[index]) ;---old ;---old ;---- for the known bodies ;---old str50 = strtrim(string(tot50, format='(E+9.2)'),2) ;---old num50 = strmid(str50,0,5) ;---old exp50 = strmid(str50,6,3) ;---old ;---old str75 = strtrim(string(tot75-tot50, format='(E+9.2)'),2) ;---old num75 = strmid(str75,0,5) ;---old exp75 = strmid(str75,6,3) ;---old fac75 = float(exp75)-float(exp50) ;---old ;---old str25 = strtrim(string(tot25-tot50, format='(E+9.2)'),2) ;---old num25 = strmid(str25,0,5) ;---old exp25 = strmid(str25,6,3) ;---old fac25 = float(exp25)-float(exp50) ;---old ;---old print, 'Total', total(stat[index,0,0]), $ ;---old float(num50), float(num25 * 10^fac25), float(num75 * 10^fac75), float(exp50), $ ;---old fracForAll[2], fracForAll[2]-fracForAll[3], fracForAll[1]-fracForAll[2], $ ;---old numN, expN, $ ;---old format='(A-25," & ",I6," & ",'+$ ;---old 'F5.2,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ . 10$^{",I2,"}$ & ",'+$ ;---old 'F6.1,"$_{\scriptscriptstyle",F+6.1,"}^{\scriptscriptstyle",F+6.1,"}$ & ",'+$ ;---old 'F3.1," . 10$^{",I1,"}$ \\")' 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