FUNCTION jpl_interp, BUF,T,NCF,NCM,NA,IFL,PV
    @compile_opt.pro        ; On error, return to caller
    common  JPL_INTERP_SAVE, PC, VC, TWOT, NP, NV

; GET CORRECT SUB-INTERVAL NUMBER FOR THIS SET OF COEFFICIENTS AND THEN GET
; NORMALIZED CHEBYSHEV TIME WITHIN THAT SUBINTERVAL.

DNA = double(NA)
DT1 = double(long(T[0]))
TC  = DNA*T[0]
L   = long(TC-DT1)

TC = 2.D0*((TC mod 1.D0)+DT1)-1.D0      ; NORMALIZED CHEBYSHEV TIME (-1<=TC<=1)

; CHECK TO SEE WHETHER CHEBYSHEV TIME HAS CHANGED, AND COMPUTE NEW POLYNOMIAL VALUES IF IT HAS.
; (PC[1] IS THE VALUE OF T1(TC) AND HENCE CONTAINS THE VALUE OF TC ON THE PREVIOUS CALL.)

NEW = n_elements(PC) EQ 0
IF NOT NEW then NEW = TC NE PC[1]

IF NEW THEN BEGIN
    TWOT  = TC+TC
    PC = dblarr(18)
    VC = dblarr(18)
    PC[0:1] = [1.0D0,TC]
    VC[0:2] = [0.0D0,1.0D0,TWOT+TWOT]       ; VC[0] is never used
    NP = 2
    NV = 3
ENDIF

PV = dblarr(NCM*IFL)                ; Output array (initialized to zero)
                                    ; At least NCF polynomials are needed
FOR I=NP,NCF-1 DO PC[I] = TWOT*PC[I-1]-PC[I-2]
NP >= NCF

; The two commented statements are basically correct, I think, but if they
; are used the test program jpl_test indicates minor differences. The do-loops
; reproduce the test cases in jpl_test exactly.

FOR I=0,NCM-1 DO BEGIN              ; Interpolate to get position for each component
    J = (L*NCM+I)*NCF
    ;PV[I] = total(PC[0:NCF-1]*BUF[J:J+NCF-1])
    FOR K=NCF-1,0,-1 DO PV[I] += PC[K]*BUF[J+K]
ENDFOR

IF IFL GE 2 THEN BEGIN
                                    ; At least NCF polynomials are needed
    FOR I=NV,NCF-1 DO VC[I] = TWOT*VC[I-1]+PC[I-1]+PC[I-1]-VC[I-2]
    NV >= NCF

    FOR I=0,NCM-1 DO BEGIN          ; Interpolate to get velocity for each component
        J = (L*NCM+I)*NCF
        ;PV[NCM+I] = total(VC[1:NCF-1]*BUF[J+1:J+NCF-1])
        FOR K=NCF-1,1,-1 DO PV[NCM+I] += VC[K]*BUF[J+K]
    ENDFOR
    PV[NCM:*] *= (DNA+DNA)/T[1]
ENDIF

RETURN, PV  &  END