; docformat = 'rst' ; ; NAME: ; readMPCORB ; PURPOSE: ; Read the MPCORB database (Minor Planet Center) and store it into an IDL structure. ; ;+ ; :Description: ; Read the MPCORB database (Minor Planet Center) and store it into an IDL structure. ; ; :Categories: ; Database, Asteroid, Orbits ; ; :Params: ; src: in, optional, type=string ; The path to the MPCORB file or an array of MPCORB lines ; ; :Returns: A structure with the orbital elements and the absolute ; magnitude. Fields are:: ; .NUM: Asteroid IAU Number ; .NAME: Asteroid Name or Provisional Designation ; .H: Absolute magnitude (Bowell HG System) ; .arcObs: Orbital arc spanned by observations ; .nbObs: Number of observations used in the computation ; .n: Mean anomaly, deg. ; .o: Argument of perihelion, deg (J2000.0) ; .Om: Longitude of ascending node, deg (J2000.0) ; .A: Semi-major axis AU ; .E: Eccentricity ; .I: Inclination (deg.) ; ; :Keywords: ; config: 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] ; mpcorb = PATH_TO_YOUR_MPCORB_FILE ; ; :Examples: ; Plot the (a,e) distribution of asteroids ; IDL> mpcorb=readMPCORB() ; IDL> plot, mpcorb.a, mpcorb.e, psym=3 ; ; :Uses: ; initIDL, epRead, date_conv ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Written in May 2014, B. Carry (IMCCE) ;- function readMPCORB, src, config=config ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Set IDL Compilation Options -------------------------------------------------- COMPILE_OPT hidden ;--I.2-- Set MPCORB File -------------------------------------------------------------- ;--I.2.1-- From default configuration if not keyword_set(src) then begin if keyword_set(config) then config= initIDL(config, /CATALOG) $ else begin confAstroIM= initIDL() confCatalog= initIDL(confAstroIM.soft.catalog, /CATALOG) endelse src = confCatalog.mpc.mpcorb ;--I.2.2-- From command-line input endif else src = strtrim(src,2) ;--I.3-- Interpret Input Source ------------------------------------------------------- srcDim=size(src) ;--I.3.1-- Orbits were Provided in Input if srcDim(0) ne 0 then begin orbits = strtrim(src,2) nbLine = srcDim[1] ;--I.3.2-- Orbits are in the Input File endif else begin ;--I.3.2/A-- Open Input File if not file_test(src,/read) then begin message, /IOERROR, 'MPCORB file not found: '+src return, -1 endif openR, idIn, src, /GET_LUN ;--I.3.2/B-- Parse Input File orbits=' ' line='' nbLine = -1L while ~EOF(idIn) do begin readf, idIn, line orbits=[orbits,strtrim(line,2)] nbLine++ endwhile orbits=orbits[1:nbLine+1] ;--I.3.2/C-- Close Input File close, idIn free_lun, idIn endelse ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Parse the Orbits from MPCORB File -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Find Valid Orbit Line -------------------------------------------------------- nbOrb = nbLine valid = where( ~strcmp(strmid(orbits,0,6),'Object') and $ ~strcmp(strmid(orbits,0,6),'') and $ ~strcmp(strmid(orbits,0,6),''), newN) if newN ne 0 then begin orbits=orbits[valid] nbOrb=newN endif ;--II.2-- Output Structure ------------------------------------------------------------- empty={pack:'', num:0L, des:'', $ H:0., G:0., $ ep:{pack:'', JD:0.D, iso:''}, $ M: 0., $ ;-Mean Anomaly peri:0., $ ;- Om:0., $ ;-Longitude of Ascending Node (Omega) n: 0., $ ;-Mean motion a: 0.D,$ ;-Semi-major axis (au) e: 0., $ ;-Eccentricity i: 0., $ ;-Inclination (degree) unc:'', $ ;-Uncertainty code ref:'', $ ;-Reference obs: {N:0L, $ ;-Number of observation first:0, $ ;-First year last:0}, $ ;-Last year opp:0L,$ ;-Number of opposition arc: 0, $ ;-Arc length (days) -- only single-opposition rms: 0. } ;-RMS residuals (") ;--II.2-- Read Each Orbit -------------------------------------------------------------- ex=replicate(empty,nbOrb) fmt='(A7,1x,F5.2,1x,F5.2,1x,A5,1x,'+$ 'F9.5,2x,F9.5,2x,F9.5,2x,F9.5,2x,F9.7,1x,F11.8,1x)' for kOrb=0, nbOrb-1 do begin ;--II.2.1-- Parse Current Orbit ex[kOrb].pack = strmid(orbits[kOrb], 0, 7 ) ;-Field 1 ex[kOrb].H = float( strmid(orbits[kOrb], 8, 5 )) ;-Field 2 ex[kOrb].G = float( strmid(orbits[kOrb], 14, 5 )) ;-Field 3 ex[kOrb].ep.pack = strmid(orbits[kOrb], 20, 5 ) ;-Field 4 ex[kOrb].M = float( strmid(orbits[kOrb], 26, 9 )) ;-Field 5 ex[kOrb].peri = float( strmid(orbits[kOrb], 37, 9 )) ;-Field 6 ex[kOrb].Om = float( strmid(orbits[kOrb], 48, 9 )) ;-Field 7 ex[kOrb].i = float( strmid(orbits[kOrb], 59, 9 )) ;-Field 8 ex[kOrb].e = float( strmid(orbits[kOrb], 70, 9 )) ;-Field 9 ex[kOrb].n = float( strmid(orbits[kOrb], 80, 11 )) ;-Field 10 ex[kOrb].a = double(strmid(orbits[kOrb], 92, 11 )) ;-Field 11 ex[kOrb].unc = strmid(orbits[kOrb], 105, 1 ) ;-Field 12 ex[kOrb].ref = strmid(orbits[kOrb], 107, 9 ) ;-Field 13 ex[kOrb].obs.N = round(float(strmid(orbits[kOrb], 117, 5 ))) ;-Field 14 ex[kOrb].opp = round(float(strmid(orbits[kOrb], 123, 3 ))) ;-Field 15 rms = strmid(orbits[kOrb], 137, 4 ) ;-Field 18 if ~strcmp(rms,' ') then ex[kOrb].rms = float(rms) try = strmid(orbits[kOrb], 132, 4 ) ;-Field 17 if strcmp(try,'days') then begin ex[kOrb].arc = round(float(strmid(orbits[kOrb], 127, 4 ))) endif else begin ex[kOrb].obs.first = round(float(strmid(orbits[kOrb], 127, 5 ))) ;-Field 16 ex[kOrb].obs.last = round(float(try)) ;-Field 17 endelse ;--II.2.2-- Convert Packed Field ex[kOrb].ep.jd = epRead( ex[kOrb].ep.pack) ex[kOrb].ep.iso = date_conv( ex[kOrb].ep.jd ) endfor return, ex end