; docformat = 'rst' ; ; NAME: ; ssoHtoV_build ; PURPOSE: ; Build the distribution of orbital elements requiered for ssoHtoV ;+ ; :Description: ; Build the distribution of orbital elements requiered for ssoHtoV ; ; :Categories: ; Database, Asteroid ; ; :Params: ; population: in, optional, type=string ; If set, only compute the dsitribution for the specified ; population (any type/subtype from http://vo.imcce.fr/webservices/skybot/?documentation) ; ; :Keywords: ; step: in, optional, type=float, default=0.25 ; The step of heliocentric distance for probability distribution ; dump: in, optional, type=string, default=(see eproc keyword) ; The directory in which the results should be stored ; plot: in, optional, type=int, default=0 ; A boolean to trigger probability plots ; eproc: in, optional, type=string, default=(from initIDL) ; The EPROC directory in the astrodata tree ; MPC: in, optional, type=string ; The path to a local version of IMCCE Discovery Information (see ssoDiscovery_build) ; configuration: in, optional, type=string ; Path to the configuration file for catalogs, at the ; minimum this file should contain the 2 following lines to be used ; by current routine:: ; [Asteroid Orbits] ; disco.imcce = PATH_TO_YOUR_IMCCE_DISCO_FILE ; eproc = PATH_TO_YOUR_BASIC_EPROC_INFORMATION ; ; :Examples: ; Build heliocentric distances for all population with a step of 0.5 AU:: ; IDL> ssoHtoV_build, step=0.5 ; ; :Uses: ; initIDL, readcol, forprint ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Written in March 2017, B. Carry (OCA) ;- pro ssoHtoV_build, population, step=step, dump=dump, plot=plot, $ MPC=MPC, eproc=eproc, configuration=configuration COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Set MPC Discovery File --------------------------------------------------------------; if not keyword_set(MPC) then begin if not keyword_set(configuration) then config = initIDL(/Catalog) $ else config = initIDL(configuration, /Catalog) MPC = config.disco.imcce endif ;--I.2-- Check Presence of MPC Discovery File ------------------------------------------------; if not file_test(MPC,/Read) then begin message, 'The MPC Discovery file cannot be read: '+strTrim(MPC,2) return endif ;--I.3-- Check Presence of Population Description File ---------------------------------------; codePATH = routine_filepath( 'ssohtov_build' ) popCSV = strMid( codePATH, 0, strlen(codePATH)-3)+'csv' ;--I.4-- Check Output Directory --------------------------------------------------------------; if not keyword_set(dump) then begin if not keyword_set(eproc) then begin if not keyword_set(configuration) then config = initIDL(/Catalog) $ else config = initIDL(configuration, /Catalog) dump = config.eproc endif else dump=eproc dump += 'Heliodist/' endif ;--I.5-- Read Population 2-D Histogram Parameters --------------------------------------------; ;--I.5.1-- Population (a,e) Range and Step readcol, popCSV, delimiter=',', format='(A,F,F,F,F,F,F)', /Silent, $ name, minA, maxA, stepA, minE, maxE, stepE, maxSun nbPop = n_elements(name) ;--I.5.2-- Fill a Structure empty={name:'',a:{min:0.,max:0.,step:0.}, e:{min:0.,max:0.,step:0.}, sun:{max:0.}} popInfo = replicate(empty,nbPop) popInfo.name = strTrim(name,2) popInfo.a.min = minA popInfo.a.max = maxA popInfo.a.step = stepA popInfo.e.min = minE popInfo.e.max = maxE popInfo.e.step = stepE popInfo.sun.max= maxSun ;--I.6-- Selection Current Population --------------------------------------------------------; if keyword_set( population ) then begin sel = where( strCmp( popInfo.name, strTrim(population,2), /Fold ) ) sel = sel[0] if sel[0] eq -1 then begin message, /ioError, 'Invalid SSO population: '+strTrim(population,2) message, /ioError, 'Population must be in: '+name return endif endif ;--I.7-- Heliocentric Step -------------------------------------------------------------------; if not keyword_set(step) then step=0.25 sufCSV='.csv' sufEPS='.eps' ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Read the MPC Discovery File -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Read the Entire File ---------------------------------------------------------------; if not keyword_set( population ) then begin disco = ssoDiscovery_read(discovery=MPC) ;--II.2-- Only Read Requested Population if *NOT* Comet ---------------------------------------; endif else if not strCmp(strLowCase(population),'comet') then begin ;--II.2.1-- Use grep to select lines - Main Population tempo='/tmp/tempo_disco.csv' spawn, 'grep "'+string(population,format='(A-9)')+'" '+MPC+' > '+tempo, cmdO, cmdE ;--II.2.2-- Check Validity - Main Population spawn, 'wc -l '+tempo, cmdO, cmdE split=strSplit( cmdO, ' ',/Extract) ;--II.2.3-- Use grep to select lines - Sub Population if strCmp(split[0], '0') then $ spawn, 'grep "'+string(population,format='(A-20)')+'" '+MPC+' > '+tempo, cmdO, cmdE ;--II.2.4-- Check Validity - Sub Population - Last Chance spawn, 'wc -l '+tempo, cmdO, cmdE split=strSplit( cmdO, ' ',/Extract) if strCmp(split[0], '0') then begin message, 'No result for the requested population: '+strTrim(population,2) spawn, 'rm '+tempo return endif ;--II.2.5-- Read the Selected Lines disco = ssoDiscovery_read(discovery=tempo) spawn, 'rm '+tempo ;--II.2.6-- Trim the popInfo Structure nbPop = 1 popInfo=popInfo[sel] endif ;--II.3-- Special Case of Comets -------------------------------------------------------------; if keyword_set(population) then if strCmp(strLowCase(population),'comet') then begin nbPop = 1 popInfo=popInfo[sel] endif ;--II.4-- Loop Over Populations --------------------------------------------------------------; for kPop=0, nbPop-1 do begin ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Build the 2-D Semi-major Axis vs Eccentricity Distribution --------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Create the 2-D Histogram --------------------------------------------------------; ;--III.1.1-- Specific Population if keyword_set( population ) then begin ;--III.1.1/A-- Comets from a Specific File if strCmp(strLowCase(population),'comet') then begin comet = readMPComet() per = where( comet.a gt 0 and finite(comet.a), nbSel ) ssoA = comet[per].a ssoE = comet[per].e print, population, n_elements(ssoA), format='(A-18,2x,I6)' ;--III.1.1/B-- Other Populations from Discovery File endif else begin ssoA = disco.a ssoE = disco.e print, population, n_elements(ssoA), format='(A-18,2x,I6)' endelse ;--III.1.2-- Looping over Populations endif else begin ;--III.1.2/A-- Comets from a Specific File if strCmp(popInfo[kPop].name,'comet',/Fold) then begin comet = readMPComet() per = where( comet.a gt 0 and finite(comet.a), nbSel ) ssoA = comet[per].a ssoE = comet[per].e print, popInfo[kPop].name, nbSel, format='(A-18,2x,I6)' ;--III.1.2/B-- Other Populations from Discovery File endif else begin case popInfo[kPop].name of 'ICB': sel = where( strCmp(disco.orbit.sub,'Inner Classical Belt',/fold), nbSel ) 'MCB': sel = where( strCmp(disco.orbit.sub,'Main Classical Belt',/fold), nbSel ) 'OCB': sel = where( strCmp(disco.orbit.sub,'Outer Classical Belt',/fold), nbSel ) 'Resonant': sel = where( strCmp(disco.orbit.sub,'Resonant',8,/fold), nbSel ) else: sel = where( strCmp(disco.orbit.main,popInfo[kPop].name,/fold) or $ strCmp(disco.orbit.sub ,popInfo[kPop].name,/fold), nbSel ) endcase print, popInfo[kPop].name, nbSel, format='(A-18,2x,I6)' ssoA = disco[sel].a ssoE = disco[sel].e endelse endelse ;--III.2-- Create the 2-D Histogram --------------------------------------------------------; h2d = hist_2d( ssoA, ssoE, $ bin1=popInfo[kPop].a.step,bin2=popInfo[kPop].e.step,$ min1=popInfo[kPop].a.min, min2=popInfo[kPop].e.min, $ max1=popInfo[kPop].a.max, max2=popInfo[kPop].e.max ) dims=size(h2d) ;--III.3-- Build the Heliocentric Distance Array -------------------------------------------; trash = ssoDistProb( 1., 0.5, delta=dSun, step=step, max=popInfo[kPop].sun.max ) probSun = dSun*0. ;--III.4-- Build the PDF of Heliocentric Distances -----------------------------------------; vecA = popInfo[kPop].a.min + findgen(dims[1]) * popInfo[kPop].a.step vecE = popInfo[kPop].e.min + findgen(dims[2]) * popInfo[kPop].e.step for kA=0, dims[1]-1 do $ for kE=0, dims[2]-1 do begin cur = ssoDistProb( vecA[kA], vecE[kE], step=step, max=popInfo[kPop].sun.max ) probSun += cur * h2d[kA,kE] endfor probSun[n_elements(probSun)-1]=0 probSun /= total(probSun) cumSun = total(probSun,/Cumulative) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Export the Distribution to a File -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Root Filename --------------------------------------------------------------------; root = popInfo[kPop].name+'-'+string(step,format='(F05.3)') ;--IV.2-- Export PDF -----------------------------------------------------------------------; forprint, dSun, probSun, cumSun, textout=dump+root+sufCSV, /Silent, $ comment='Heliodist, Probability, Cumulative', $ format='(F8.3,",",E10.4,",",E10.4)' ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Plot the Distribution -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(plot) then begin ;--V.1-- Open Graphic Environment --------------------------------------------------------; cgPS_open, Filename=dump+root+sufEPS, /metric, /decomposed, /portrait, $ xsize=30, ysize=20, language_level=2, /quiet ;--V.2-- Inner Solar System Plot ---------------------------------------------------------; if popInfo[kPop].a.max le 10 then begin ;--V.2.1-- PDF xRange=[0,10] cgPlot, dSun, probSun, xRange=xRange, psym=10, $ xStyle=1, xTitle='Heliocentric distance (AU)', $ yStyle=8, yTitle='Probability (/'+strTrim(string(step,format='(F8.3)'),2)+' AU)' ;--V.2.2-- Cumulative PDF cgPlot, dSun, cumSun, xRange=xRange, /NoErase, $ color='Cornflower Blue', $ xStyle=5, yStyle=5 ;--V.3-- Outer Solar System Plot ---------------------------------------------------------; endif else begin ;--V.3.1-- PDF xRange = [1,1000] cgPlot, dSun, probSun, xRange=xRange, psym=10, /xLog, $ xStyle=1, xTitle='Heliocentric distance (AU)', $ yStyle=8, yTitle='Probability (/'+strTrim(string(step,format='(F8.3)'),2)+' AU)' ;--V.3.2-- Cumulative PDF cgPlot, dSun, cumSun, xRange=xRange, /xLog, /NoErase, $ color='Cornflower Blue', $ xStyle=5, yStyle=5 endelse ;--V.4-- Plot the Axis for PDF of Heliocentric Distance ----------------------------------; cgAxis, yAxis=1, yStyle=1, yRange=[0,1], yTitle='Cumulative probability', color='Cornflower Blue' ;--V.5-- Open Graphic Environment --------------------------------------------------------; cgPS_close, /png, Delete_PS=0 endif endfor end