sbtrp - dec. 1975 - r. goldstein
velocities from numerical differentiation added  (nvel= -1)
by j.f.chandler, 1977 may
this is a control routine for interpolating a body from its
integration tape (or n-body tape)
input parameters:
jd        integer value of the coord. time of the desired
point
fract     fractional part of the coordinate time of the point
nvel  = 0 positions determined
1 velocities determined
-1 velocities determined
.gt.1 treated as 1
lcntl = 0 use results of last call if possible
1 do not use results (must recalculate y vectors)
-1 same as 0, except that nb1 is forced to 2 (this may
require some otherwise unnecessary y-vector setups)
lcntl must be set to 1 the first time the routine is called
output is put in xsb vector

enter here to get acceleration for same (jd,fract) as last
call. for everett method, uses whichever y vectors are
available and corresponding differentiating interpolation.
for hermite method, will have to calculate new y vector.