;; ; ; Compute the completeness of spectral/photometric surveys VS MPC ; ; astorb = readAstorb() ; save=astorb ; w={a:astorb.a, h:astorb.h} ; astorb=w binH = 0.5 rangH = [2,22] hMin = fltarr(3) hMax = fltarr(3) nbBin = (rangH(1)-rangH(0))/binH + 1 rangA = [1.780, 2.501] hMin(0) = 14.0 & hMax(0) = 16 ; hMin(0) = 12. & hMax(0) = 14.5 pv=0.178885 rho=2000.*1e9 samp1 = where( astorb.a ge rangA(0) and astorb.a lt rangA(1), nbSam1 ) histoH = histogram( astorb.h(samp1), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) cumulH = total( histoH, /CUMUL ) mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) hSet = where( hArr ge hMin(0) and hArr le hMax(0) ) lin = linfit( hArr(hSet), alog10(cumulH(hSet)) ) modelMB = lin(0) + hArr * lin(1) cumulN = 10.^modelMB totalMB=modelMB*0. for kB=nbBin-1,1,-1 do $ totalMB(kB) = cumulN(kB) - cumulN(kB-1) lowH=where(hArr le hMax(0)) fracMB = histoH/totalMB fracMB(lowH) = 1 modelMB=totalMB modelMB(lowH) = histoH[lowH] cgplot, hArr, cumulH,/ylog,yr=[1,1e6] cgPlot, hArr, cumulN, /over, color='grey' cgplot, hArr, histoH cgPlot, hArr, totalMB, /over, color='grey' cgPlot, hArr, fracMB,psym=10 stop ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ; compile_opt idl2 cwd = '/home/bcarry/data/mining/sdss/' readDisco = 0 readMOC = 0 readDMC = 1 readLarge = 1 ;--II.1-- Database paths --------------------------------------------------------------------- conf= initIDL() confCata = initIDL(conf.soft.catalog, /catalog) pathDisco = confCata.disco.imcce pathMOC = cwd+'ADR4.dat' pathDMC = cwd+'alluniq.dat' pathLarge = cwd+'inputs/totalclean.txt' ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Read Large Databases -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(verbose) then print, '--- Reading databases' ;--II.1-- MPC Discovery Information ---------------------------------------------------------- if readDisco eq 1 then begin disco=ssoDiscovery_read(discovery=pathDisco, /dump) if keyword_set(verbose) then print, ' ... MPC discovery read' endif ;--II.2-- SDSS Moving Object Catalog --------------------------------------------------------- if readMOC eq 1 then begin ;--II.2.1-- Read the whole MOC moc=sdssReadMOC(pathMOC, release=4) ;--II.2.2-- Select only once each object with number/provisional designation sel = where( moc.ssoNum ne 0 or ~strcmp(moc.ssoName,'-',1), nbSel ) numName=string(moc[sel].ssoNum,format='(I06)')+strtrim(moc[sel].ssoName,2) uniqSel=uniq(numName,sort(numName)) moc=moc[sel[uniqSel]] if keyword_set(verbose) then print, ' ... SDSS MOC read' endif ;--II.3-- DeMeo-Carry SDSS Taxonomy ---------------------------------------------------------- if readDeMeoCarry eq 1 then begin sdss = sdssReadTaxonomy( pathDMC ) if keyword_set(verbose) then print, ' ... DeMeo-Carry Taxonomy read' endif ;--II.4-- Remaining Large SSO from Spectral Surveys ----------------------------------------- if readLarge eq 1 then begin readfmt, pathLarge, 'I8,A20,F9.2,F8.3,A6,1x,A10,F8.4,1x,A10', $ num, name, h, a, type, class, alb, source, /silent large = { num:num, des:name, h: h, a:a, cVal:type, cSrc:class, qVal:alb, qSrc:source } if keyword_set(verbose) then print, ' ... Totalclean read' endif ;--II.5-- Read AstOrb ----------------------------------------------------------------------- if readLarge eq 1 then begin readfmt, pathLarge, 'I8,A20,F9.2,F8.3,A6,1x,A10,F8.4,1x,A10', $ num, name, h, a, type, class, alb, source, /silent large = { num:num, des:name, h: h, a:a, cVal:type, cSrc:class, qVal:alb, qSrc:source } if keyword_set(verbose) then print, ' ... Totalclean read' endif ;-----------old ; readJDD = 0 ; readPDS = 0 ; readMOC = 0 ; readUNIQ= 0 ; readASTORB= 0 ; ; searchMIS = 0 ; window, 0 ; window, 2 ; root = '/home/bcarry/data/mining/sdss/' ; if readJDD eq 1 then begin ; readfmt, '/astrodata/mpc/discovery.tab', $ ; '(I6,1x,A20,1x,F5.2,1x,D9.1,1x,F9.6,1x,F10.8,1x,F12.8)', $ ; num, des, h, jdd, a, e, i, /SILENT ; ; mpINFO = { num: num, $ ; des: des, $ ; h: h, $ ; jdd: jdd, $ ; a: a, $ ; e: e, $ ; i: i } ; ; print, ' -- Discovery read --' ; endif ; j1981 = date_conv( '1981-01-01T00:00:00', 'JULIAN' ) ; j1982 = date_conv( '1982-01-01T00:00:00', 'JULIAN' ) ; j1983 = date_conv( '1983-01-01T00:00:00', 'JULIAN' ) ; ; v1981 = where( mpINFO.jdd le j1981, nb1981) ; v1982 = where( mpINFO.jdd le j1982, nb1982) ; v1983 = where( mpINFO.jdd le j1983, nb1983) ; ; print, nb1981, nb1982, nb1983 ; ; if readPDS eq 1 then begin ; readfmt, root+'cescafile/totalclean.txt', $ ; 'I8,A20,F9.2,F8.3,A6,1x,A10,F8.4,1x,A10', $ ; num, name, h, a, type, class, alb, source, $ ; /silent ; pdsINFO = { num: num, des: name, $ ; h: h, a:a, $ ; cVal: type, cSrc: class, $ ; qVal: alb, qSrc: source } ; ; print, ' -- Totalclean read --' ; endif ; if readMOC eq 1 then begin ; allmoc=read_moc(root+'cescafile/ADR4.dat', 1) ; nbMOC = float(n_elements( allmoc.moID )) ; ; mocDYNA = return_struct(allmoc,'moc') ; object number/orbit ; mocDYNA.name = strtrim(mocDYNA.name,2) ; ; mocNum = where( mocDYNA.num ne 0, nbN ) ; mocPrv = where( mocDYNA.name ne '-' and mocDYNA.num eq 0, nbP ) ; ; uniqMOCn = uniq( mocDYNA.num(mocNum), sort(mocDYNA.num(mocNum)) ) ; uniqMOCp = uniq( mocDYNA.name(mocPrv), sort(mocDYNA.name(mocPrv)) ) ; ; mocINFO = { num: mocDYNA.num([mocNum(uniqMOCn), mocPrv(uniqMOCp)]), $ ; des: mocDYNA.name([mocNum(uniqMOCn), mocPrv(uniqMOCp)]), $ ; h: mocDYNA.hmag([mocNum(uniqMOCn), mocPrv(uniqMOCp)]), $ ; a: mocDYNA.a([mocNum(uniqMOCn), mocPrv(uniqMOCp)]), $ ; e: mocDYNA.e([mocNum(uniqMOCn), mocPrv(uniqMOCp)]), $ ; i: mocDYNA.i([mocNum(uniqMOCn), mocPrv(uniqMOCp)]) } ; ; ;; for k=0, n_elements( pdsINFO.num )-1 do begin ;; pp=where( pdsINFO.num(k) eq mocINFO.num ) ;; if pp(0) ne -1 then mocINFO.h(pp(0)) = 99 ;; endfor ; ; ; print, ' -- SDSS MOC ADR4 read --' ; endif rangA = [1.780, 2.501] binH = 0.5 rangH = [2,22] ;help, allmoc,/str samp1 = where( allmoc.a_sdss ge rangA(0) and allmoc.a_sdss lt rangA(1), nbSam1 ) histoH = histogram( allmoc.hmag(samp1), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) cumulH = total( histoH, /CUMUL ) plot, hArr, cumulH, psym=10, /ylog, yr=[1,1e6] stop if readUNIQ eq 1 then begin uniqINFO = read_uniq( root+'alluniq.dat' ) for k=0, n_elements( pdsINFO.num )-1 do begin pp=where( pdsINFO.num(k) eq uniqINFO.num ) if pp(0) ne -1 then uniqINFO.h(pp(0)) = 99 endfor print, ' -- UNIQ read --' endif if readASTORB eq 1 then begin astOrbFile = '/astrodata/ASTORB/astorb.dat' astorb =read_astorb( astOrbFile ) ;-Update PDS with MPC H numb=where( pdsINFO.num ne 0 ) pdsINFO.h(numb) = astorb.h(pdsINFO.num(numb)-1) ;-Update MOC with MPC H numb=where( mocINFO.num ne 0 ) mocINFO.h(numb) = astorb.h(mocINFO.num(numb)-1) print, ' -- AstOrb read --' endif rho=2000.*1e9 hMin = fltarr(3) hMax = fltarr(3) binH = 0.5 rangH = [2,22] nbBin = (rangH(1)-rangH(0))/binH + 1 histoMB = lonarr( nbBin, 3) cumulMB = lonarr( nbBin, 3) modelMB = fltarr( nbBin, 3) totalMB = fltarr( nbBin, 3) massMB = fltarr( nbBin, 3) massBus = fltarr( nbBin, 3) massMOC = fltarr( nbBin, 3) massUniq = fltarr( nbBin, 3) surveys = fltarr( nbBin, 3) sdssMOC = fltarr( nbBin, 3) uniqMOC = fltarr( nbBin, 3) fitArr = fltarr(2,3) h1 = 8.0 h2 = 10.0 ; diam = [5,1.,0.5,0.1] ; nbDiam = n_elements( diam) ; inputDiam = fltarr( nbDiam,3 ) ; numberDiam = fltarr( nbDiam,3 ) ;--------------------------------------------------------------- ;-IMB rangA = [1.780, 2.501] hMin(0) = 14.0 & hMax(0) = 16 ; hMin(0) = 12. & hMax(0) = 14.5 pv=0.178885 samp1 = where( mpINFO.a ge rangA(0) and mpINFO.a lt rangA(1), nbSam1 ) histoH = histogram( mpINFO.h(samp1), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) cumulH = total( histoH, /CUMUL ) ; mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * total( histoH * 10.^(-3*hArr/5.) ) mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) histoMB(*,0) = histoH cumulMB(*,0) = cumulH massMB(*,0) = mass hSet = where( hArr ge hMin(0) and hArr le hMax(0) ) lin = linfit( hArr(hSet), alog10(cumulH(hSet)) ) modelMB(*,0) = lin(0) + hArr * lin(1) fitArr(*,0) = [ lin(0),lin(1)] ; hDiam = -5.*alog10( diam*sqrt(pv)/1329. ) ; inputDiam(*,0)= hDiam ; numberDiam(*,0) = lin(0) + hDiam * lin(1) samp2 = where( pdsINFO.a ge rangA(0) and pdsINFO.a lt rangA(1), nbSam2 ) histoH = histogram( pdsINFO.h(samp2), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) surveys(*,0) = histoH mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) massBus(*,0) = mass samp3 = where( mocINFO.a ge rangA(0) and mocINFO.a lt rangA(1), nbSam3 ) histoH = histogram( mocINFO.h(samp3), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) sdssMOC(*,0) = histoH mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) massMOC(*,0) = mass samp4 = where( uniqINFO.a ge rangA(0) and uniqINFO.a lt rangA(1), nbSam4 ) histoH = histogram( uniqINFO.h(samp4), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) uniqMOC(*,0) = histoH mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) massUniq(*,0) = mass if searchMIS eq 1 then begin pp = where( mpINFO.h(samp1) ge h1 and mpINFO.h(samp1) lt h2) list1 = mpINFO.des(samp1(pp)) nums1 = mpINFO.num(samp1(pp)) hmag1 = mpINFO.h(samp1(pp)) pp = where( pdsINFO.h(samp2) ge h1 and pdsINFO.h(samp2) lt h2) list2 = pdsINFO.des(samp2(pp)) n1 = n_elements( list1 ) n2 = n_elements( list2 ) print, 'IMB - missing: ',n1-n2 for kk=0, n1-1 do begin test = where( list1(kk) eq list2 ) if test(0) eq -1 then print, nums1(kk), list1(kk), hmag1(kk) endfor endif ;--------------------------------------------------------------- ;-MMB rangA = [2.501, 2.825] hMin(1) = 13 & hMax(1) = 15 ; hMin(1) = 14.0 & hMax(1) = 16 pv=0.141372 samp1 = where( mpINFO.a ge rangA(0) and mpINFO.a lt rangA(1), nbSam1 ) histoH = histogram( mpINFO.h(samp1), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) cumulH = total( histoH, /CUMUL ) ; mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * total( histoH * 10.^(-0.6*hArr) ) mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-0.6*hArr) ) histoMB(*,1) = histoH cumulMB(*,1) = cumulH massMB(*,1) = mass hSet = where( hArr ge hMin(1) and hArr le hMax(1) ) lin = linfit( hArr(hSet), alog10(cumulH(hSet)) ) modelMB(*,1) = lin(0) + hArr * lin(1) fitArr(*,1) = [ lin(0),lin(1)] ; hDiam = -5.*alog10( diam*sqrt(pv)/1329. ) ; inputDiam(*,1)= hDiam ; numberDiam(*,1) = lin(0) + hDiam * lin(1) samp2 = where( pdsINFO.a ge rangA(0) and pdsINFO.a lt rangA(1), nbSam2 ) histoH = histogram( pdsINFO.h(samp2), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) surveys(*,1) = histoH mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) massBus(*,1) = mass samp3 = where( mocINFO.a ge rangA(0) and mocINFO.a lt rangA(1), nbSam3 ) histoH = histogram( mocINFO.h(samp3), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) sdssMOC(*,1) = histoH mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) massMOC(*,1) = mass samp4 = where( uniqINFO.a ge rangA(0) and uniqINFO.a lt rangA(1), nbSam4 ) histoH = histogram( uniqINFO.h(samp4), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) uniqMOC(*,1) = histoH mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) massUniq(*,1) = mass if searchMIS eq 1 then begin pp = where( mpINFO.h(samp1) ge h1 and mpINFO.h(samp1) lt h2) list1 = mpINFO.des(samp1(pp)) nums1 = mpINFO.num(samp1(pp)) hmag1 = mpINFO.h(samp1(pp)) pp = where( pdsINFO.h(samp2) ge h1 and pdsINFO.h(samp2) lt h2) list2 = pdsINFO.des(samp2(pp)) n1 = n_elements( list1 ) n2 = n_elements( list2 ) print, 'MMB - missing: ',n1-n2 for kk=0, n1-1 do begin test = where( list1(kk) eq list2 ) if test(0) eq -1 then print, nums1(kk), list1(kk), hmag1(kk) endfor endif ;--------------------------------------------------------------- ;-OMB rangA = [2.825, 3.278] hMin(2) = 12. & hMax(2) = 14.5 ; hMin(2) = 14.0 & hMax(2) = 16 pv=0.0903660 & tit='OMB' samp1 = where( mpINFO.a ge rangA(0) and mpINFO.a lt rangA(1), nbSam1 ) histoH = histogram( mpINFO.h(samp1), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) cumulH = total( histoH, /CUMUL ) ; mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * total( histoH * 10.^(-0.6*hArr) ) mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-0.6*hArr) ) histoMB(*,2) = histoH cumulMB(*,2) = cumulH massMB(*,2) = mass hSet = where( hArr ge hMin(2) and hArr le hMax(2) ) lin = linfit( hArr(hSet), alog10(cumulH(hSet)) ) modelMB(*,2) = lin(0) + hArr * lin(1) fitArr(*,2) = [ lin(0),lin(1)] ; hDiam = -5.*alog10( diam*sqrt(pv)/1329. ) ; inputDiam(*,2)= hDiam ; numberDiam(*,2) = lin(0) + hDiam * lin(1) samp2 = where( pdsINFO.a ge rangA(0) and pdsINFO.a lt rangA(1), nbSam2 ) histoH = histogram( pdsINFO.h(samp2), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) surveys(*,2) = histoH mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) massBus(*,2) = mass samp3 = where( mocINFO.a ge rangA(0) and mocINFO.a lt rangA(1), nbSam3 ) histoH = histogram( mocINFO.h(samp3), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) sdssMOC(*,2) = histoH mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) massMOc(*,2) = mass samp4 = where( uniqINFO.a ge rangA(0) and uniqINFO.a lt rangA(1), nbSam4 ) histoH = histogram( uniqINFO.h(samp4), bins=binH, min=rangH(0), max=rangH(1), loc=hArr ) uniqMOC(*,2) = histoH mass = (!PI*rho/6.) * ((1329.D/sqrt(pv))^3.) * ( histoH * 10.^(-3*hArr/5.) ) massUniq(*,2) = mass if searchMIS eq 1 then begin pp = where( mpINFO.h(samp1) ge h1 and mpINFO.h(samp1) lt h2) list1 = mpINFO.des(samp1(pp)) nums1 = mpINFO.num(samp1(pp)) hmag1 = mpINFO.h(samp1(pp)) pp = where( pdsINFO.h(samp2) ge h1 and pdsINFO.h(samp2) lt h2) list2 = pdsINFO.des(samp2(pp)) n1 = n_elements( list1 ) n2 = n_elements( list2 ) print, 'OMB - missing: ',n1-n2 for kk=0, n1-1 do begin test = where( list1(kk) eq list2 ) if test(0) eq -1 then print, nums1(kk), list1(kk), hmag1(kk) endfor endif ;--------------------------------------------------------------- ;-PLOTS device, decom=0 loadct, 39, /SILENT plot, hArr, cumulMB(*,0), xr=[10,22], yr=[2e2,3e5],/YLOG oplot, hArr, cumulMB(*,1) oplot, hArr, cumulMB(*,2) oplot, hArr, 10.^modelMB(*,0), color=230 oplot, hArr, 10.^modelMB(*,1), color=230 oplot, hArr, 10.^modelMB(*,2), color=230 diff = modelMB ; alog10( cumulMB ) cumulN = 10.^diff for kB=nbBin-1,1,-1 do begin totalMB(kB,0) = cumulN(kB,0) - cumulN(kB-1,0) totalMB(kB,1) = cumulN(kB,1) - cumulN(kB-1,1) totalMB(kB,2) = cumulN(kB,2) - cumulN(kB-1,2) endfor fracMB = histoMB/totalMB fracMB(where(hArr le hMax(0)),0) = 1 fracMB(where(hArr le hMax(1)),1) = 1 fracMB(where(hArr le hMax(2)),2) = 1 MBR = 4 if MBR eq 0 then begin print, ' Inner main belt' wset, 0 plot, hArr, histoMB(*,0), xr=[3,22], yr=[1.1,3e5],/YLOG, xst=1,yst=1,psym=10 oplot, hArr, totalMB(*,0), color=230 oplot, hArr, surveys(*,0), color=130, psym=10 oplot, hArr, sdssMOC(*,0), color=90, psym=10 oplot, hArr, uniqMOC(*,0), color=90, linestyl=1, psym=10 wset, 2 plot, hArr, fracMB(*,0), xr=[3,22], yr=[0,1.1], xst=1,yst=1 oplot, hArr, fracMB(*,0)*surveys(*,0)/histoMB(*,0), psym=10, color=230 oplot, hArr, fracMB(*,0)*sdssMOC(*,0)/histoMB(*,0), psym=10, color=80 oplot, hArr, fracMB(*,0)*uniqMOC(*,0)/histoMB(*,0), psym=10, color=90, linestyl=1 print, total(uniqMOC(*,0)), total(sdssMOC(*,0)), total(uniqMOC(*,0))/float(total(sdssMOC(*,0))) combine = fracMB(*,0)*(surveys(*,0)+sdssMOC(*,0))/histoMB(*,0) combine( where(combine ge 1)) =1 ; oplot, hArr, combine, psym=10, color=130, linestyle=2 combine = fracMB(*,0)*(surveys(*,0)+uniqMOC(*,0))/histoMB(*,0) combine( where(combine ge 1)) =1 oplot, hArr, combine, psym=10, color=130, linestyle=1 endif if MBR eq 1 then begin print, ' Middle main belt' wset, 0 plot, hArr, histoMB(*,1), xr=[3,22], yr=[1.1,3e5],/YLOG, xst=1,yst=1, psym=10 oplot, hArr, totalMB(*,1), color=230 oplot, hArr, surveys(*,1), color=130, psym=10 oplot, hArr, sdssMOC(*,1), color=90, psym=10 oplot, hArr, uniqMOC(*,1), color=90, linestyl=1, psym=10 wset, 2 plot, hArr, fracMB(*,1), xr=[3,22], yr=[0,1.1], xst=1,yst=1 oplot, hArr, fracMB(*,1)*surveys(*,1)/histoMB(*,1), psym=10, color=230 oplot, hArr, fracMB(*,1)*sdssMOC(*,1)/histoMB(*,1), psym=10, color=80 oplot, hArr, fracMB(*,1)*uniqMOC(*,1)/histoMB(*,1), psym=10, color=90, linestyl=1 print, total(uniqMOC(*,1)), total(sdssMOC(*,1)), total(uniqMOC(*,1))/float(total(sdssMOC(*,1))) combine = fracMB(*,1)*(surveys(*,1)+sdssMOC(*,1))/histoMB(*,1) combine( where(combine ge 1)) =1 ; oplot, hArr, combine, psym=10, color=130, linestyle=2 combine = fracMB(*,1)*(surveys(*,1)+uniqMOC(*,1))/histoMB(*,1) combine( where(combine ge 1)) =1 oplot, hArr, combine, psym=10, color=130, linestyle=1 endif if MBR eq 2 then begin print, ' Outer main belt' wset, 0 plot, hArr, histoMB(*,2), xr=[3,22], yr=[1.1,3e5],/YLOG, xst=1,yst=1 oplot, hArr, totalMB(*,2), color=230 oplot, hArr, surveys(*,2), color=130 oplot, hArr, sdssMOC(*,2), color=90 oplot, hArr, uniqMOC(*,2), color=90, linestyl=1 wset, 2 plot, hArr, fracMB(*,2), xr=[3,22], yr=[0,1.1], xst=1,yst=1 oplot, hArr, fracMB(*,2)*surveys(*,2)/histoMB(*,2), psym=10, color=230 oplot, hArr, fracMB(*,2)*sdssMOC(*,2)/histoMB(*,2), psym=10, color=80 oplot, hArr, fracMB(*,2)*uniqMOC(*,2)/histoMB(*,2), psym=10, color=90, linestyl=1 print, total(uniqMOC(*,2)), total(sdssMOC(*,2)), total(uniqMOC(*,2))/float(total(sdssMOC(*,2))) combine = fracMB(*,2)*(surveys(*,2)+sdssMOC(*,2))/histoMB(*,2) combine( where(combine ge 1)) =1 ; oplot, hArr, combine, psym=10, color=130, linestyle=2 combine = fracMB(*,2)*(surveys(*,2)+uniqMOC(*,2))/histoMB(*,2) combine( where(combine ge 1)) =1 oplot, hArr, combine, psym=10, color=130, linestyle=1 endif ; window, 3 ; plot, hArr, ( 1329/sqrt(0.2) ) * 10^(-hArr/5.), /YLOG, $ ; xr=[3,22], yr=[0.1,1000], xst=1,yst=1 ; h=[5,10,15,20.] forprint, h, ( 1329/sqrt(0.2) ) * 10^(-h/5.), textout=2 h=[16,15,14.5] forprint, h, ( 1329/sqrt(0.2) ) * 10^(-h/5.), textout=2 d=[0.1,1,10,100.] forprint, d, -5*alog10( d*sqrt(0.2)/1329.), textout=2 wdelete set_plot, 'PS' device, xo=0, yo=0 device, xsize=19.4, ysize=9.0 device, /ENCAPSULATED device, /COLOR device, BITS_PER_PIXEL=8 device, filename=root+'sdss-completeness.eps' !P.thick=3 xBot=1200 yBot=900 xLen=6000 yLe1=2500 yLe2=5000 xSep= 50 ySep= 50 xRang = [5,21.5] yRangH= [1.1,3e5] yRangF= [0,1.05] !P.charsize=0.8 !P.charthick=1.5 yTickH = '10!U'+string(indgen(5)+1,format='(I1)')+'!D' yTickF = string(indgen(6)*20,format='(I3)') yTitH = 'Distribution' yTitF = 'Fraction (%)' xTit = ' ' Title = ['Inner', 'Middle', 'Outer'] cInd = [ 0, 130, 230, 80, 80] cTab = [ 0, 0, 39, 39, 39] cSty = [ 0, 2, 0, 0, 1] hExtra = where( hArr ge 11 ) for kM=0, 2 do begin ; print, kM if kM eq 1 then xTit= 'Absolute magnitude (H)' if kM ge 1 then yTickH=replicate(' ',10) loadct, 39, /SILENT plot, hArr, histoMB(*,kM), xr=xRang, yr=yRangH, /YLOG, xst=9,yst=1, $ position= [xBot+ kM*(xLen+xSep), $ yBot+yLe1+ySep, $ xBot+ (kM+1)*(xLen)+ kM*xSep, $ yBot+yLe1+ySep+yLe2 ], /DEVICE, /NOERASE, /NODATA, $ xthick=2, ythick=2, $ xtickname=replicate(' ',4), $ ytickname=yTickH, $ ytitle=yTitH, color=cInd(0), linestyl=cSty(0), $ xticklen =0.025 axis, xaxis=1, xr=(1329/sqrt(0.2))*10^(-xRang/5.), /xlog, xst=1, xthick=2 oplot, hArr, histoMB(*,kM), color=cInd(0), linestyl=cSty(0) oplot, hArr(hExtra), totalMB(hExtra,kM), color=cInd(1), linestyl=cSty(1) oplot, hArr, surveys(*,kM), color=cInd(2), linestyl=cSty(2) oplot, hArr, sdssMOC(*,kM), color=cInd(3), linestyl=cSty(3) oplot, hArr, uniqMOC(*,kM), color=cInd(4), linestyl=cSty(4) ; oplot, hArr, 10^(0.3*hArr), thick=4 xyouts, total(xRang)/2., 1e5, /DATA, Title(kM), align=0.5 if kM eq 0 then begin xKey = 5.7 kLen = 1.5 yKey = 1e3 * exp(findgen(5)) oplot, xKey+[0,kLen], yKey(0)+[0,0], color=cInd(4), linestyl=cSty(4) oplot, xKey+[0,kLen], yKey(1)+[0,0], color=cInd(3), linestyl=cSty(3) oplot, xKey+[0,kLen], yKey(2)+[0,0], color=cInd(2), linestyl=cSty(2) oplot, xKey+[0,kLen], yKey(3)+[0,0], color=cInd(1), linestyl=cSty(1) oplot, xKey+[0,kLen], yKey(4)+[0,0], color=cInd(0), linestyl=cSty(0) yKey = 8e2 * exp(findgen(5)) xyouts, xKey+1.3*kLen, yKey(4), 'MPC' xyouts, xKey+1.3*kLen, yKey(3), 'Power law' xyouts, xKey+1.3*kLen, yKey(2), 'Spectra' xyouts, xKey+1.3*kLen, yKey(1), 'SDSS MOC' xyouts, xKey+1.3*kLen, yKey(0), 'SDSS cut' endif plot, hArr, fracMB(*,kM), xr=xRang, yr=yRangF, xst=1,yst=1, $ position= [xBot+ kM*(xLen+xSep), $ yBot, $ xBot+ (kM+1)*(xLen)+ kM*xSep, $ yBot+yLe1 ], /DEVICE, /NOERASE, /NODATA, $ xthick=2, ythick=2, $ ytickname=yTickF, ytitle=yTitF, xtitle=xTit, color=cInd(0), linestyl=cSty(0), $ xticklen =0.025*yLe2/yLe1 ;tag oplot, hArr, loadct, 39, /SILENT oplot, xRang, [1,1], color=cInd(1), linestyl=cSty(1), thick=2 loadct, 0, /SILENT final = fracMB(*,kM)*(surveys(*,kM)+uniqMOC(*,kM))/histoMB(*,kM) pF = where( hArr ge xRang(0)+.5 and hArr le xRang(1) and final ge yRangF(0), nbPM ) loadct, 0, /silent for kF=0, nbPM-1 do begin polyfill, hArr(pF(kF))+binH*[-1,1,1,-1]/2., $ [final(pF(kF)),final(pF(kF)),0,0], color=230 endfor polyfill, xRang(0)+[0,binH, binH, 0], $ [1,1,0,0], color=230 oplot, hArr, final, psym=10, color=200, linestyl=cSty(0), thick=2 loadct, 39, /silent oplot, hArr, fracMB(*,kM), psym=10, color=cInd(0), linestyl=cSty(0), thick=2 oplot, hArr, fracMB(*,kM)*surveys(*,kM)/histoMB(*,kM), psym=10, color=cInd(2), linestyl=cSty(2), thick=2 oplot, hArr, fracMB(*,kM)*surveys(*,kM)/histoMB(*,kM), psym=10, color=cInd(2), linestyl=cSty(2), thick=2 oplot, hArr, fracMB(*,kM)*sdssMOC(*,kM)/histoMB(*,kM), psym=10, color=cInd(3), linestyl=cSty(3), thick=2 oplot, hArr, fracMB(*,kM)*uniqMOC(*,kM)/histoMB(*,kM), psym=10, color=cInd(4), linestyl=cSty(4), thick=2 if kM eq 0 then begin loadct, 0, /silent xKey = 17 yKey = .70 xKLen = 1.25 yKLen = .20 polyfill, xKey+[0,1,1,0]*xKLen, yKey+[1,1,0,0]*yKLen, color=230 oplot, xKey+[0,1,1,0,0]*xKLen, yKey+[1,1,0,0,1]*yKLen, color=200 xyouts, xKey+1.3*xKLen, yKey+.35* yKLen, 'Total' endif plot, hArr, fracMB(*,kM), xr=xRang, yr=yRangF, xst=1,yst=1, $ position= [xBot+ kM*(xLen+xSep), $ yBot, $ xBot+ (kM+1)*(xLen)+ kM*xSep, $ yBot+yLe1 ], /DEVICE, /NOERASE, /NODATA, $ ytickname=yTickF, ytitle=yTitF, xtitle=xTit, color=cInd(0), linestyl=cSty(0) ; print, total(uniqMOC(*,kM)), total(sdssMOC(*,kM)), total(uniqMOC(*,kM))/float(total(sdssMOC(*,kM))) ; xyouts, xBot + kM*(xLen+xSep) + 0.5*xLen, yBot+yLe1+yLe2+4*ySep, /DEVICE, $ ; Title(kM), align=0.5 yTickF=replicate(' ',6) yTitH=' ' yTitF=' ' xTit = ' ' ; endfor xyouts, 250, yBot+yLe1+yLe2+4*ySep, /DEVICE, 'Diameter:' xyouts, 1*xLen+ 50, yBot+yLe1+yLe2+4*ySep, /DEVICE, '(km)' device, /CLOSE set_plot, 'X' forprint, title, hMin, hMax, total(massMB,1), -5*fitArr(1,*), textout=2, $ format='(A-8,F5.2,1x,F5.2,4x,E8.2,4x,F5.2)' ; histoMB = lonarr( nbBin, 3) ; cumulMB = lonarr( nbBin, 3) ; modelMB = fltarr( nbBin, 3) ; totalMB = fltarr( nbBin, 3) ; massMB = fltarr( nbBin, 3) ; ; surveys = fltarr( nbBin, 3) ; sdssMOC = fltarr( nbBin, 3) ; uniqMOC = fltarr( nbBin, 3) print, total(massMB) print, total(massBus) print, total(massMOC) print, total(massUniq) print, 100*[ total(massBus), total(massMOC), total(massUniq)]/ total(massMB) massMB = total( massMB,2 ) relatM = massMB/total(massMB) cumulM = total(massMB,/CUMUL)/total(massMB) print, '' print, 'mass per bin' forprint, hArr, relatM, cumulM, textout=2 print, '' print, 'number of objects' forprint, hArr, $ totalMB(*,0), $ totalMB(*,1), $ totalMB(*,2), total(totalMB,2), textout=2, subset=where(hArr ge 17) ; for k =0,2 do begin ; print, k+1 ; print, inputDiam(*,k) ; print, 10.^numberDiam(*,k) ; endfor end