doELT = 0 ;-Search orbital elements readAO= 01 ;-Read astorb doSRC = 0 ;-Run Greenstreet source mapper doCAT = 01 ;-Concatenate all mapping results doSUM = 0 ;-Combine NEOs/Class into Src/Class probabilities nS=5 srcLabel=['nu6', '3:1', 'M-C', 'OMB', 'JFC'] classType = ['A', 'B', 'C', 'D', 'K', 'L', 'Q', 'S', 'V', 'X'] nC = n_elements(classType) root = '/home/bcarry/data/mining/mithneos/' dirSRC = root+'dynamics/source_region/' dirELT = dirSRC+'elements/' dirMAP = dirSRC+'mapper/' srcMapper = root+'greenstreet/Release_NoIntegrations/restime/separated/source_mapper/source_mapper_ex' readcol, root+'inputRPB/MITHNEOS_2015_08_14_table.csv', num, name, comet, des, $ orb,dv, $ visData,nirData, visRef,visRefSym,nirRef,nirRefSym, $ taxRef,otherRef,wave, dataRef,t1,t2,class, $ format='(L,A,A,A,A,F,'+$ 'A,A,A,A,A,A,'+$ 'A,A,A,A,A,A,A)', $ delim=',', skipline=1, /silent nbAst = n_elements(des) ; forprint, num, name, des, class, textout=2, format='(I6,1x,A-12,1x,A-12,1x,A-6)', subset=indgen(10) ;--------------------------------------------------------------------------------; ;--------------------------------------------------------------------------------; ;-- I -- Retrieve Asteroid Orbital Elements ----------------------------; ;--------------------------------------------------------------------------------; ;--------------------------------------------------------------------------------; if doELT eq 1 then begin if readAO eq 1 then astorb=readAstorb() for kAst=0, nbAst-1 do begin ;-name or designation if strcmp(name(kAst),'') or strcmp(strtrim(name(kAst),2),'0') then cur=des(kAst) else cur=name(kAst) Id=string(num(kAst),format='(I06)')+'-'+strtrim(string(cur,format='(A-22)'),2)+'.dat' ;-clean spaces if strmid(id,11,1) eq ' ' then strput, id, '_', 11 space = strpos(id,' ') if space ne -1 then strput, id, '_', space weird = strpos(id,"'") if weird ne -1 then strput, id, '_', weird ;-retrieve orbital elements if num(kAst) ne 0 then begin aei = astorb(num(kAst)-1) endif else begin p = where( strcmp( cur, astorb.name,/fold ) ) if p(0) ne -1 then aei = astorb(p(0)) else begin print, 'Issue: '+cur endelse endelse forprint, strtrim(string([aei.a, aei.e, aei.i], format='(F)'),2), format='(A-12)', $ /nocomment, textout=dirELT+Id, /silent endfor endif ;--------------------------------------------------------------------------------; ;--------------------------------------------------------------------------------; ;-- II -- Run Greenstreet Source Region Model ----------------------------; ;--------------------------------------------------------------------------------; ;--------------------------------------------------------------------------------; if doSRC eq 1 then begin for kAst=0, nbAst-1 do begin ;-name or designation if strcmp(name(kAst),'') or strcmp(strtrim(name(kAst),2),'0') then cur=des(kAst) else cur=name(kAst) Id=string(num(kAst),format='(I06)')+'-'+strtrim(string(cur,format='(A-22)'),2)+'.dat' ;-clean spaces if strmid(id,11,1) eq ' ' then strput, id, '_', 11 space = strpos(id,' ') if space ne -1 then strput, id, '_', space weird = strpos(id,"'") if weird ne -1 then strput, id, '_', weird ;-User Greenstreet src mapper spawn, 'cd '+file_dirname(srcmapper)+' && ./'+$ file_basename(srcMapper)+' < '+dirELT+Id spawn, 'mv '+file_dirname(srcmapper)+'/source_mapper.txt '+dirMap+Id endfor endif ;--------------------------------------------------------------------------------; ;--------------------------------------------------------------------------------; ;-- III -- Concatenation of Source Region Results -------------------------; ;--------------------------------------------------------------------------------; ;--------------------------------------------------------------------------------; if doCAT eq 1 then begin empty = {nu6:0., m31:0., MC:0., OB:0., JFC:0.} src=replicate( {val:empty, unc:empty}, nbAst ) for kAst=0, nbAst-1 do begin ;-name or designation if strcmp(name(kAst),'') or strcmp(strtrim(name(kAst),2),'0') then begin name(kAst) = '' cur=des(kAst) endif else cur=name(kAst) Id=string(num(kAst),format='(I06)')+'-'+strtrim(string(cur,format='(A-22)'),2)+'.dat' ;-clean spaces if strmid(id,11,1) eq ' ' then strput, id, '_', 11 space = strpos(id,' ') if space ne -1 then strput, id, '_', space weird = strpos(id,"'") if weird ne -1 then strput, id, '_', weird ;-Read source mapping results readfmt, dirMap+Id, '15x,F5.3,4x,F5.3',p,dp,/silent src(kAst).val.nu6 = p(27) src(kAst).val.m31 = p(28) src(kAst).val.mc = p(29) src(kAst).val.ob = p(30) src(kAst).val.jfc = p(31) src(kAst).unc.nu6 = dp(27) src(kAst).unc.m31 = dp(28) src(kAst).unc.mc = dp(29) src(kAst).unc.ob = dp(30) src(kAst).unc.jfc = dp(31) endfor forprint, num, des, name, class, $ src.val.nu6, src.val.m31, src.val.MC, src.val.OB, src.val.JFC, $ src.unc.nu6, src.unc.m31, src.unc.MC, src.unc.OB, src.unc.JFC, $ format='(I6,1x,A-12,1x,A-12,1x,A-8,1x,5(F6.4,1x),1x,5(F6.4,1x))', $ textout=dirMAP+'summary.dat', /SILENT, /NOCOMMENT endif ;--------------------------------------------------------------------------------; ;--------------------------------------------------------------------------------; ;-- IV -- Class vs Source Region -------------------------; ;--------------------------------------------------------------------------------; ;--------------------------------------------------------------------------------; ;---- if doSUM eq 1 then begin ;---- ;---- ;-READ source mapping summary file ;---- readfmt, dirMAP+'summary.dat', $ ;---- 'I6,1x,A-12,1x,A-3,1x,F6.4,1x,F6.4,1x,F6.4,1x,F6.4,1x,F6.4,1x,1x,F6.4,1x,F6.4,1x,F6.4,1x,F6.4,1x,F6.4', $ ;---- sumNum, sumDesig, neaClass, $ ;---- Vnu6, Vm31, VMC, VOB, VJFC, $ ;---- Unu6, Um31, UMC, UOB, UJFC, $ ;---- /silent ;---- neaClass=strtrim(neaClass,2) ;---- ;---- empty = {nu6:0., m31:0., MC:0., OB:0., JFC:0.} ;---- src=replicate( {val:empty, unc:empty}, nbAst ) ;---- src.val.nu6 = Vnu6 ;---- src.val.m31 = Vm31 ;---- src.val.mc = Vmc ;---- src.val.ob = Vob ;---- src.val.jfc = Vjfc ;---- ;---- src.unc.nu6 = Unu6 ;---- src.unc.m31 = Um31 ;---- src.unc.mc = Umc ;---- src.unc.ob = Uob ;---- src.unc.jfc = Ujfc ;---- ;---- ;---- ;-Normalize src probablilites ;---- vec = [[src.val.nu6], [src.val.m31], [src.val.mc], [src.val.ob], [src.val.jfc]] ;---- unc = [[src.unc.nu6], [src.unc.m31], [src.unc.mc], [src.unc.ob], [src.unc.jfc]] ;---- ;---- norm=vec ;---- noru=unc ;---- ;---- tot = total(vec,1) ;---- megaTot = total(vec) ;---- ;---- ;---- print, '# source & total', tot, total(tot) ;---- print, '% source ', 100.*tot/total(tot) ;---- ;---- ;---- for k=0,nS-1 do begin ;---- norm(*,k) = vec(*,k)/tot(k) ;---- noru(*,k) = unc(*,k)/tot(k) ;---- endfor ;---- ;---- ;---- ;---- ;---- ;---- ;-Count fraction of Class per source region ;---- clasSrc = fltarr(nC,nS,2) ;---- srcNorm = fltarr(nC,nS,2) ;---- normZ = fltarr(nS) ;---- for kS=0,nS-1 do begin ;---- for kC=0,nC-1 do begin ;---- ;---- cur = where( strcmp(neaClass, classType(kC),/fold), nbCur ) ;---- ;---- tot1 = total( vec(cur,kS) ) ;---- tot2 = nbCur*total( vec(cur,kS)/unc(cur,kS), /NAN ) / total( 1./unc(cur,kS), /NAN ) ;---- tot3 = nbCur*meanWithUnc( vec(cur,kS), unc(cur,kS) ) ;---- ;---- print, srcLabel(kS), class(kC), strtrim(string(nbCur),2), $ ;---- tot1,tot2, tot3(0),format='(A-4,2x,A1,2x,I4,4x,3(F7.2,2x))' ;---- ;---- std1 = sqrt( total( unc(cur,kS)^2 ) ) ;---- ;---- tot=tot3(0) ;---- std=tot3(1) ;---- ;---- tot=tot1 ;---- std=std1 ;---- ;---- ;---- clasSrc(kC,kS,*) = tot ;---- clasSrc(kC,kS,1) = std ;---- ;---- endfor ;---- ;---- normZ(kS) = 1.*total( clasSrc(*,kS,0), /NAN ) ;---- ;---- for kC=0,nC-1 do begin ;---- srcNorm(kC,kS,0) = clasSrc(kC,kS,0) / normZ(kS) ;---- srcNorm(kC,kS,1) = clasSrc(kC,kS,1) / normZ(kS) ;---- ;---- cur = where( strcmp(neaClass, classType(kC),/fold), nbCur ) ;---- if nbCur eq 0 then begin ;---- srcNorm(kC,kS,*) = 0 ;---- endif else begin ;---- ;---- ssMarg1 = 1.96 * sqrt( srcNorm(kC,kS,0) * (1-srcNorm(kC,kS,0)) /nbCur ) ;-95% confidence ;---- ssMarg2 = 0.98 * sqrt( srcNorm(kC,kS,0) * (1-srcNorm(kC,kS,0)) /nbCur ) ;-95% confidence ;---- ssMarg3 = 0.98 / sqrt( nbCur ) ;-95% confidence ;---- ;---- FPC = sqrt( (megaTot - nbCur) / (megaTot-1) ) ;---- ;----; print,nbCur, srcNorm(kC,kS,0), ssMarg1, ssMarg2, ssMarg3, FPC, 0.98/sqrt(megaTot) ;---- ;---- ssMarg = ssMarg2 * FPC + 0.98/sqrt(megaTot) ;---- ;----; print, 100*[srcNorm(kC,kS,0), srcNorm(kC,kS,1), ssMarg] ;---- srcNorm(kC,kS,1) = sqrt( srcNorm(kC,kS,1)^2 + ssMarg^2 ) ;---- ;---- endelse ;---- ;---- endfor ;---- ;---- endfor ;---- ;---- endif ;- ;- ;- ;- ;- dirMap = root+'ancillary/src/' ;- dirIn = dirMap ;- dirEx = dirMap+'outputs/' ;- ;- ;- CASE set of ;- 'DSS': begin ;-;-obspm nea = read_uniq(root+'class/svo-dss-noz-uniq.dat') ;- nea = read_uniq(root+'class/svo-dss-noz-uniq.dat') ;- end ;- ELSE: ;- ENDCASE ;- nbAst = n_elements(nea.num) ;- ;- ;- ;- nC=10 ;- nS=5 ;- ;- class = ['A ', 'B ', 'C ', 'D ', 'K ', 'L ', 'Q ', 'S ', 'V ', 'X '] ;- ;- ;- ;- empty = {nu6:0., m31:0., MC:0., OB:0., JFC:0.} ;- ;- ;- src=replicate( {val:empty, unc:empty}, nbAst ) ;- ;- ;- if doIN eq 1 then $ ;- for k=0, nbAst-1 do begin ;- Id=string(nea.num(k),format='(I06)')+'-'+strtrim(string(nea.desig(k),format='(A-12)'),2)+'.dat' ;- if strmid(id,11,1) eq ' ' then strput, id, '_', 11 ;- ;- space = strpos(id,' ') ;- if space ne -1 then strput, id, '_', space ;- ;- forprint, strtrim(string([nea.a(k), nea.e(k), nea.i(k)], format='(F)'),2), format='(A-12)', /nocomment, textout=dirIN+Id, /silent ;- endfor ;- ;- ;- ;- if doCAT eq 1 then begin ;- for k=0, nbAst-1 do begin ;- Id=string(nea.num(k),format='(I06)')+'-'+strtrim(string(nea.desig(k),format='(A-12)'),2)+'.dat' ;- if strmid(id,11,1) eq ' ' then strput, id, '_', 11 ;- ;- space = strpos(id,' ') ;- if space ne -1 then strput, id, '_', space ;- ;- readfmt, dirEX+id, '15x,F5.3,4x,F5.3',p,dp,/silent ;- ;- src(k).val.nu6 = p(27) ;- src(k).val.m31 = p(28) ;- src(k).val.mc = p(29) ;- src(k).val.ob = p(30) ;- src(k).val.jfc = p(31) ;- ;- src(k).unc.nu6 = dp(27) ;- src(k).unc.m31 = dp(28) ;- src(k).unc.mc = dp(29) ;- src(k).unc.ob = dp(30) ;- src(k).unc.jfc = dp(31) ;- ;- endfor ;- ;- forprint, nea.num, nea.desig, nea.class, $ ;- src.val.nu6, src.val.m31, src.val.MC, src.val.OB, src.val.JFC, $ ;- src.unc.nu6, src.unc.m31, src.unc.MC, src.unc.OB, src.unc.JFC, $ ;- format='(I6,1x,A-12,1x,A-3,1x,5(F6.4,1x),1x,5(F6.4,1x))', $ ;- textout=dirEX+'summary.dat', /SILENT, /NOCOMMENT ;- endif ;- ;- ;- ;- ;- ;- ;- ;- if doSR eq 1 then begin ;- ;- readfmt, dirEX+'summary.dat', $ ;- 'I6,1x,A-12,1x,A-3,1x,F6.4,1x,F6.4,1x,F6.4,1x,F6.4,1x,F6.4,1x,1x,F6.4,1x,F6.4,1x,F6.4,1x,F6.4,1x,F6.4', $ ;- num, desig, neaClass, $ ;- Vnu6, Vm31, VMC, VOB, VJFC, $ ;- Unu6, Um31, UMC, UOB, UJFC, $ ;- /silent ;- ;- src.val.nu6 = Vnu6 ;- src.val.m31 = Vm31 ;- src.val.mc = Vmc ;- src.val.ob = Vob ;- src.val.jfc = Vjfc ;- ;- src.unc.nu6 = Unu6 ;- src.unc.m31 = Um31 ;- src.unc.mc = Umc ;- src.unc.ob = Uob ;- src.unc.jfc = Ujfc ;- ;- vec = [[src.val.nu6], [src.val.m31], [src.val.mc], [src.val.ob], [src.val.jfc]] ;- unc = [[src.unc.nu6], [src.unc.m31], [src.unc.mc], [src.unc.ob], [src.unc.jfc]] ;- ;-; unc(where(unc eq 0)) = 0.01 ;- ;- ;- norm=vec ;- noru=unc ;- ;- tot = total(vec,1) ;- megaTot = total(vec) ;- ;- ;- print, '# source & total', tot, total(tot) ;- print, '% source ', 100.*tot/total(tot) ;- ;- ;- for k=0,nS-1 do begin ;- norm(*,k) = vec(*,k)/tot(k) ;- noru(*,k) = unc(*,k)/tot(k) ;- endfor ;- ;- ;- srcLabel=['nu6', '3:1', 'M-C', 'OMB', 'JFC'] ;- ;- ;- clasSrc = fltarr(nC,nS,2) ;- srcNorm = fltarr(nC,nS,2) ;- normZ = fltarr(nS) ;- for kS=0,nS-1 do begin ;- ;- for kC=0,nC-1 do begin ;- ;- cur = where( neaClass eq class(kC), nbCur ) ;- ;- tot1 = total( vec(cur,kS) ) ;- tot2 = nbCur*total( vec(cur,kS)/unc(cur,kS), /NAN ) / total( 1./unc(cur,kS), /NAN ) ;- tot3 = nbCur*meanWithUnc( vec(cur,kS), unc(cur,kS) ) ;- ;- print, srcLabel(kS), class(kC), strtrim(string(nbCur),2), $ ;- tot1,tot2, tot3(0),format='(A-4,2x,A1,2x,I4,4x,3(F7.2,2x))' ;- ;- std1 = sqrt( total( unc(cur,kS)^2 ) ) ;- ;- ;- tot=tot3(0) ;- std=tot3(1) ;- ;- tot=tot1 ;- std=std1 ;- ;- ;- clasSrc(kC,kS,*) = tot ;- clasSrc(kC,kS,1) = std ;- ;- endfor ;- ;- normZ(kS) = 1.*total( clasSrc(*,kS,0), /NAN ) ;- ;- for kC=0,nC-1 do begin ;- srcNorm(kC,kS,0) = clasSrc(kC,kS,0) / normZ(kS) ;- srcNorm(kC,kS,1) = clasSrc(kC,kS,1) / normZ(kS) ;- ;- cur = where( neaClass eq class(kC),nbCur ) ;- if nbCur eq 0 then begin ;- srcNorm(kC,kS,*) = 0 ;- endif else begin ;- ;- ssMarg1 = 1.96 * sqrt( srcNorm(kC,kS,0) * (1-srcNorm(kC,kS,0)) /nbCur ) ;-95% confidence ;- ssMarg2 = 0.98 * sqrt( srcNorm(kC,kS,0) * (1-srcNorm(kC,kS,0)) /nbCur ) ;-95% confidence ;- ssMarg3 = 0.98 / sqrt( nbCur ) ;-95% confidence ;- ;- FPC = sqrt( (megaTot - nbCur) / (megaTot-1) ) ;- ;-; print,nbCur, srcNorm(kC,kS,0), ssMarg1, ssMarg2, ssMarg3, FPC, 0.98/sqrt(megaTot) ;- ;- ssMarg = ssMarg2 * FPC + 0.98/sqrt(megaTot) ;- ;-; print, 100*[srcNorm(kC,kS,0), srcNorm(kC,kS,1), ssMarg] ;- srcNorm(kC,kS,1) = sqrt( srcNorm(kC,kS,1)^2 + ssMarg^2 ) ;- ;- endelse ;- ;- endfor ;- ;- ;- ;- endfor ;- ;- ;- ;- ;- ;- ;-;-obspm readcol, '/obs/bcarry/data/mining/neas/ancillary/num35.txt', $ ;- readcol, dirMap+'num35.txt', $ ;- nId, nClass, nAll, n100, n50, n20, n5, n3, format='(A,A,F,F,F,F,F,F)', /Silent ;- g = where( nClass ne 'R' and nId ne 'i') ;- frac3 = 100. * n3(g)/total(n3(g)) ;- ;- ;- ;- ;- ;-;------------- CONVERT TO PERCENT HERE!!!! ;- ;- ;- srcNorm *= 100. ;- print, ' ', srcLabel,format='(A3,5(5x,A3,5x))' ;- forprint, class, $ ;- srcNorm(*,0,0),srcNorm(*,0,1), $ ;- srcNorm(*,1,0),srcNorm(*,1,1), $ ;- srcNorm(*,2,0),srcNorm(*,2,1), $ ;- srcNorm(*,3,0),srcNorm(*,3,1), $ ;- srcNorm(*,4,0),srcNorm(*,4,1), $ ;- n3(g), frac3, $ ;- format='(A3,1x,5(2(F5.2,1x),1x),I4,1x,F5.2)', textout=2 ;- ;- sep="' & '" ;- forprint, class, $ ;- srcNorm(*,0,0),srcNorm(*,0,1), $ ;- srcNorm(*,1,0),srcNorm(*,1,1), $ ;- srcNorm(*,2,0),srcNorm(*,2,1), $ ;- srcNorm(*,3,0),srcNorm(*,3,1), $ ;- srcNorm(*,4,0),srcNorm(*,4,1), $ ;- n3(g), frac3, $ ;- format='(A3,'+sep+',5(2(F4.1,'+sep+')),I4,'+sep+',F4.1," \\")', textout=2 ;- ;- ;- ;- ;- ;- ;- compLabel =['A','B','C','D','K','L','S','V','X'] ;- compId =[0,1,2,3,4,5,7,8,9] ;- ;- ;- ;- for SrId=0, nS-1 do begin ;-; srId = 01 ;- ;- x = srcNorm(compId,srId,0) ;- x(6) += srcNorm(6,srId,0) ;add Q to S ;- y = frac3(compId) ;- ;- a=regress(x,y, const=b, corr=correl, chisq=chi2, yfit=line) ;- print, '--------------'+srcLabel(srId) ;- print, 'coef :', a ;- print, 'origin:', b ;- print, 'correl:', correl ;- print, 'chi2 :', chi2 ;- ;- plot, x,y, psym=4, /iso ;- oplot, x, line ;- ;-;wait,1 ;- endfor ;- ;- ;-; ;-; print, total( srcNorm(*,0,0)) ;-; print, total( srcNorm(*,1,0)) ;-; print, total( srcNorm(*,2,0)) ;-; print, total( srcNorm(*,3,0)) ;-; print, total( srcNorm(*,4,0)) ;-; ;- ;-;------------------------------ figure histogram ;- ;- fIMB = frac3(compId) ;- ;- fNu6 = srcNorm(compId,0,0) ;- fNu6(6) += srcNorm(6,0,0) ;add Q to S ;- ;- f3_1 = srcNorm(compId,1,0) ;- f3_1(6) += srcNorm(6,1,0) ;add Q to S ;- ;- fMC = srcNorm(compId,2,0) ;- fMC(6) += srcNorm(6,2,0) ;add Q to S ;- ;- fOB = srcNorm(compId,3,0) ;- fOB(6) += srcNorm(6,3,0) ;add Q to S ;- ;- fJF = srcNorm(compId,4,0) ;- fJF(6) += srcNorm(6,4,0) ;add Q to S ;- ;- ;- fracs = [ [fIMB], [fNu6], [f3_1], [fMC], [fOB], [fJF] ] ;-; help, fracs ;- ;- xArr = indgen(9)+1 ;- xShi = 0.21*[-1,0,1,2] ;- bins = 0.1 ;- names = ['Inner Main Belt 3-5 km', $ ;- 'NEAs from '+cgGreek('nu',/PS)+'!D6!N', $ ;- 'NEAs from 3:1 MMR', $ ;- 'NEAs from M-C space', $ ;- 'NEAs from Outer Belt', $ ;- 'NEAs from JFC'] ;- ;- cols=['Cornflower Blue', $ ;- 'Goldenrod', $ ;- 'Almond', $ ;- 'Orange', $ ;- 'Gray', $ ;- 'Navy'] ;- ;- wdelete ;- ;- set_plot, 'PS' ;- device, xsize=18, ysize=13 ;- device, filename=dirMap+'fig-sources.eps' ;- device, /TIMES ;- device, /ENCAPSULATED ;- !P.FONT=0 ;- ;- plot, [0], [0], /NODATA, /DEVICE, $ ;- xst=1, xr=[0,10], xminor=1, $ ;- xtickint=1, xtickname=[' ',compLabel,' '], $ ;- yst=1, yr=[0,60], yminor=4, $ ;- ytickint=20, $ ;- ytitle='Fraction of the population (%)', $ ;- position=[1300,800,17500,12500] ;- ;- ;- xKey = 0.5 ;- yKey = 55 ;- kShi = 5 ;- ;- for kS=0,2 do begin ;- for kC=0,8 do begin ;- ;- xArr = (kC+1) + xShi(kS) + bins*[1,1,-1,-1] ;- yArr = fracs(kC,kS) * [0,1,1,0] ;- ;- print, compLabel(kC), fracs(kC,0)-fracs(kC,1), fracs(kC,0)-fracs(kC,2), $ ;- format='(A1,2x,F4.1,2x,F4.1)' ;- ;- polyfill, xarr, yarr, color=cgcolor(cols(kS)),/DATA ;- ;- endfor ;- ;- xArr = xKey + bins*[1,1,-1,-1]*1.5 ;- yArr = yKey + bins*[1,-1,-1,1]*1.5*8 - kShi*kS ;- ;- polyfill, xarr, yarr, color=cgcolor(cols(kS)),/DATA ;- xyOuts, xKey + 2.5*bins, yKey-kShi*kS-.8, names(kS) ;- ;- ;- endfor ;- ;- ;- ;- device, /close ;- set_plot, 'X' ;- ;- ;- endif ;- ;- ;-end ;- end