; read a Minor Planet Circular file

function mpc_rd2, nelem, body=body, fn=fn, nskip=nskip

    if not keyword_set(nskip) then begin
       nskip = 0
    endif
    if not keyword_set(fn) then begin
       fn = '~/projects/plseason09/20111101_regime/mpc_w_header.txt'
       nskip = 4
    endif

    form = '('+ $
      ' A7, 1X'   + $;  Number or provisional designation (in packed form)
      ',F5.2, 1X' + $;  Absolute magnitude, H
      ',F5.2, 1X' + $;  Slope parameter, G
      ',A5, 1X'   + $;  Epoch (in packed form, .0 TT)
      ',F9.5, 2X' + $;  Mean anomaly at the epoch, in degrees
      ',F9.5, 2X' + $;  Argument of perihelion, J2000.0 (degrees)
      ',F9.5, 2X' + $;  Longitude of the ascending node, J2000.0 (degrees)
      ',F9.5, 2X' + $;  Inclination to the ecliptic, J2000.0 (degrees)
      ',F9.7, 1X' + $;  Orbital eccentricity
      ',F11.8, 1X'+ $;  Mean daily motion (degrees per day)
      ',F11.7, 2X'+ $;  Semimajor axis (AU)
      ',A1, 1X'   + $;  Uncertainty parameter, U
      ',A9, 1X'   + $;  Reference
      ',A5, 1X'   + $;  Number of observations
      ',I3, 1X'   + $;  Number of oppositions
      ',I4    '   + $;  Year of first observation or arclen (days)
      ',A1    '   + $;  ',-'                       or empty
      ',A4, 1X'   + $;  Year of last observation  or 'days'
      ',F4.2, 1X' + $;  r.m.s residual (")
      ',A3, 1X'   + $;  Coarse indicator of perturbers
      ',A3, 1X'   + $;  Precise indicator of perturbers
      ',A10, 1X'  + $ ;  Computer name
      ',Z4, 1X'   + $ ;  4-hexdigit flags
      ',A27, 1X'  + $ ;  Readable designation
      ',A8'       + $ ;  Date of last observation included in orbit solution
      ')'

    readfmt, fn, form, $
      id, H, G, Epoch_packed, M, APF, LAN, i, e, n, A, $
      u, ref, nobs, nopp, yrfirst, dash, yrlast, rms, pertc, pertp, comp, $
      flags, desig, lastobs, /silent, skip=nskip
    nelem = n_elements(id)

    elem1 = { id:'', elem:dblarr(7), type:0, desig:'', $
              h:0., g:0., u:'', ref:'', nobs:'', nopp:0, yrfirst:0, yrlast:'', $
              rms:0.0, pertc:'', pertp:'', comp:'', flags:' ', lastobs:0L}
    elems = replicate(elem1, nelem)
    for ielem = 0, nelem-1 do begin
        elems[ielem].id = id[ielem]
        elems[ielem].h = h[ielem]
        elems[ielem].g = g[ielem]
        epoch = jdparse(unpack_mpcepoch(Epoch_packed[ielem]))
        tau = epoch - m[ielem]/n[ielem]       ; Julian Date of last periapsis
        elems[ielem].elem[0] = e[ielem]
        elems[ielem].elem[1] = i[ielem]
        elems[ielem].elem[2] = LAN[ielem]
        elems[ielem].elem[3] = APF[ielem]
        elems[ielem].elem[4] = tau
        elems[ielem].elem[5] = n[ielem]
        elems[ielem].elem[6] = a[ielem]
        elems[ielem].type = (flags[ielem]) and (64-1) ; 
;                       10  Centaur
;                       14  Plutino
;                       15  Other resonant TNO
;                       16  Cubewano
;                       17  Scattered disk

        elems[ielem].ref = ref[ielem]
        elems[ielem].nobs = nobs[ielem]
        elems[ielem].nopp = nopp[ielem]
        elems[ielem].yrfirst = yrfirst[ielem]
        elems[ielem].yrlast = yrlast[ielem]
        elems[ielem].rms = rms[ielem]
        elems[ielem].pertc = pertc[ielem]
        elems[ielem].pertp = pertp[ielem]
        elems[ielem].comp = comp[ielem]
        elems[ielem].flags = flags[ielem]
        elems[ielem].desig = desig[ielem]
        elems[ielem].lastobs = lastobs[ielem]

     end

    if n_elements(body) ne 0 then begin
       i = matchfirst(1,strmatch(desig,'*'+body+'*'))
       if i eq -1 then begin
          nelem = 0
          return, 0
       endif else begin
          nelem = 1
          return, elems[i]
       endelse
    endif

    return, elems

end