subroutine lgpv9(tnorm,xt,yt,zt,vxt,vyt,vzt,
     +    xout,yout,zout,vxout,vyout,vzout)
      implicit real*8 (a-h,o-z)
      dimension xt(9),yt(9),zt(9),vxt(9),vyt(9),vzt(9),c(9),t(9)
      data t/0.0d0,1.0d0,2.0d0,3.0d0,4.0d0,5.0d0,6.0d0,7.0d0,8.0d0/
      data tol/1.0d-10/
      data c/40320.0d0,-5040.0d0,1440.0d0,-720.0d0,576.0d0,
     +  -720.0d0,1440.0d0,-5040.0d0,40320.0d0/
c.....................................................................
c                                                                    .
c                                                                    .
c   Subroutine: lgpv9  Author: B. Remondi  Date: Jan/15/1983         .
c                                                                    .
c   Function:  lgpv9 is a 9-th order Lagrangian interpolater.        .
c              It has been adapted, here, quite severely, for        .
c              speed.  To extend this adaptation to 11-th order      .
c              would be straightforward.  The routine is inten-      .
c              tionally nongeneral for max. speed.  The indepen-     .
c              dent variable, tnorm, is normalized time.  Thus       .
c              the position and velocity arrays (e.g., xt(k))        .
c              hold values for tnorm = 0,1,2,...,8, respectively.    .
c              The equal spacing between array values has allowed    .
c              the c(j), j=1,..,9, to be computed.  This algoritm    .
c              can be made yet faster.  By the way, for tnorm        .
c              requests on (or within tol) the integer times         .
c              t=0,1,...,8 interpolation is avoided.                 .
c                                                                    .
c   Inputs: tnorm - Normalized time.  Thus the user must scale his   .
c                   time before calling this subroutine.  This       .
c                   subroutine will work even for normalized time    .
c                   outside of the range  -1.LE.tnorm.LE.9, however  .
c                   wild results are possible.                       .
c                                                                    .
c           xt(k),  yt(k),  zt(k)  -  position in any units.         .
c           vxt(k), vyt(k), vzt(k) -  velocity in any units.         .
c                                                                    .
c   Outputs: xout, yout, zout, vxout, vyout, vzout in the same units..
c                                                                    .
c                                                                    .
c.....................................................................
c                                                                    .
c
c  First see if tnorm=t(i), i=1,9.  If so, we are at an exact time
c  of the stored data, so interpolation is not needed.
c
      intgrt=tnorm+tol
c
c     write(6,121) tnorm
c 121 format(' "lgpv9" tnorm',f20.15)
c
      if(dabs(intgrt-tnorm).lt.2.0d0*tol) go to 1000
 50   prodn=1.d0
      xout=0.d0
      yout=0.d0
      zout=0.d0
      vxout=0.d0
      vyout=0.d0
      vzout=0.d0
      do 100 i=1,9
      prodn=(tnorm-t(i))*prodn
      prodd=1.0d0/(c(i)*(tnorm-t(i)))
      xout=xt(i)*prodd+xout
      yout=yt(i)*prodd+yout
      zout=zt(i)*prodd+zout
      vxout=vxt(i)*prodd+vxout
      vyout=vyt(i)*prodd+vyout
 100  vzout=vzt(i)*prodd+vzout
      xout=xout*prodn
      yout=yout*prodn
      zout=zout*prodn
      vxout=vxout*prodn
      vyout=vyout*prodn
      vzout=vzout*prodn
      return
 1000 continue
c
c To reach this code, tnorm = an integer value (within tol).
c We must never violate the inequality   0.LE.intgrt.LE.8 .
c
      itrubl=iabs(intgrt-4)
      if(itrubl.gt.4) go to 50
      xout=xt(intgrt+1)
      yout=yt(intgrt+1)
      zout=zt(intgrt+1)
      vxout=vxt(intgrt+1)
      vyout=vyt(intgrt+1)
      vzout=vzt(intgrt+1)
      return
      end