; docformat = 'rst' ; ; NAME: ; readMPCOMET ; PURPOSE: ; Read the MPCOMET database (Minor Planet Center) and store it into an IDL structure. ; ;+ ; :Description: ; Read the MPCOMET database (Minor Planet Center) and store it into an IDL structure. ; The file can be found at https://www.minorplanetcenter.net/data "CometEls.txt" ; ; :Categories: ; Database, Comet, Orbits ; ; :Params: ; src: in, optional, type=string ; The path to the MPCOMET file or an array of MPCOMET lines ; ; :Returns: A structure with the orbital elements and the absolute ; magnitude. Fields are:: ; .NUM: Comet IAU Number ; .NAME: Comet IAU 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] ; mpcomet = PATH_TO_YOUR_MPCOMET_FILE ; ; :Examples: ; Plot the (a,e) distribution of asteroids ; IDL> MPCOMET=readMPCOMET() ; IDL> cgPlot, MPCOMET.a, MPCOMET.e, psym='Filled circle',/xLog, xRange=[1,10000] ; ; :Uses: ; initIDL, date_conv ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Written in June 2017, B. Carry (OCA) ;- function readMPCOMET, src, config=config COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Set MPCOMET File -------------------------------------------------------------------- ;--I.1.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.mpcomet ;--I.1.2-- From command-line input endif else src = strTrim(src,2) ;--I.2-- Interpret Input Source ------------------------------------------------------------- srcDim=size(src) ;--I.2.1-- Orbits were Provided in Input if srcDim[0] ne 0 then begin orbits = strTrim(src,2) nbLine = srcDim[1] ;--I.2.2-- Orbits are in the Input File endif else begin ;--I.2.2/A-- Open Input File if not file_test(src,/read) then begin message, /ioError, 'MPCOMET file not found: '+src return, -1 endif openR, idIn, src, /GET_LUN ;--I.2.2/B-- Parse Input File orbits=' ' line='' nbLine = -1L while ~EoF(idIn) do begin readf, idIn, line orbits=[orbits,line] nbLine++ endwhile orbits=orbits[1:nbLine+1] ;--I.2.2/C-- Close Input File close, idIn free_lun, idIn endelse ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Parse the Orbits from MPCOMET File -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Find Valid Orbit Line -------------------------------------------------------------- nbOrb = nbLine ;--II.2-- Output Structure ------------------------------------------------------------------- empty={num:0L, name:'', des:'', type:'', $ tpp:{JD:0.D, iso:''}, $ ;-Time of passage of perihelion ep:{JD:0.D, iso:''}, $ ;-Epoch of perturbed solution a: 0.D,$ ;-Semi-major axis (au) e: 0., $ ;-Eccentricity p: 0., $ ;-Orbital period peri:0., $ ;-Perihelion argPeri:0., $ ;-Argument of perihelion node:0., $ ;-Longitude of Ascending Node (degrees) i:0., $ ;-Inclination (degrees) H:0., G:0., $ ;-Absolute magnitude and slope ref:'' } ;-Reference ;--II.3-- Read Each Orbit -------------------------------------------------------------------- ex=replicate(empty,nbOrb) for kOrb=0, nbOrb-1 do begin ;--II.2.1-- Comet Designation str = strMid(orbits[kOrb],0,4) if not strCmp(str,' ') then ex[kOrb].num = round(float(str)) ex[kOrb].type = strMid(orbits[kOrb], 4, 1) ex[kOrb].des = strTrim(strMid(orbits[kOrb], 5, 7),2) ex[kOrb].name = strTrim(strMid(orbits[kOrb],102,30),2) ;--II.2.2-- Time of Passage to Perihelion year = strMid( orbits[kOrb], 14, 4 ) month = strMid( orbits[kOrb], 19, 2 ) day = strMid( orbits[kOrb], 22, 2 ) frac = float(strMid( orbits[kOrb], 24, 5 )) hh = floor( frac*24 ) mm = floor(( frac*24 - hh )* 60) ss = (( frac*24 - hh )* 60 - mm )*60 ex[kOrb].tpp.iso = year+'-'+month+'-'+day+'T'+$ string(hh,format='(I02)')+':'+$ string(mm,format='(I02)')+':'+$ string(ss,format='(F07.4)') ex[kOrb].tpp.jd = date_conv(ex[kOrb].tpp.iso,'JULIAN') ;--II.2.3-- Orbital Elements ex[kOrb].peri = float(strMid( orbits[kOrb], 31, 9 )) ex[kOrb].e = float(strMid( orbits[kOrb], 41, 8 )) ex[kOrb].a = ex[kOrb].peri / ( 1 - ex[kOrb].e ) ex[kOrb].argPeri = float(strMid( orbits[kOrb], 51, 8 )) ex[kOrb].node = float(strMid( orbits[kOrb], 61, 8 )) ex[kOrb].i = float(strMid( orbits[kOrb], 71, 8 )) ;--II.2.4-- Epoch for perturbed solution year = strMid( orbits[kOrb], 81, 4 ) month = strMid( orbits[kOrb], 85, 2 ) day = strMid( orbits[kOrb], 87, 2 ) if strCmp(year,' ') then begin ex[kOrb].ep.iso = ex[kOrb-1].ep.iso endif else begin ex[kOrb].ep.iso = year+'-'+month+'-'+day+'T00:00:00.000' endelse ex[kOrb].ep.jd = date_conv(ex[kOrb].ep.iso,'JULIAN') ;--II.2.5-- Magnitude ex[kOrb].H = float(strMid( orbits[kOrb], 91, 4 )) ex[kOrb].G = float(strMid( orbits[kOrb], 96, 4 )) ;--II.2.6-- Reference ex[kOrb].ref = strTrim(strMid( orbits[kOrb], 159, 9 ),2) ; print, kOrb ; print, ex[kOrb].ep.iso ; print, ex[kOrb].tpp.iso endfor return, ex end