C+
C NAME:
C	jpl_interp
C PURPOSE:
C	Differentiates and interpolates a set of Chebyshev
C	coefficients to give position and velocity
C CALLING SEQUENCE:
	subroutine jpl_interp(BUF,T,NCF,NCM,NA,IFL,pv)
C INPUTS:
C	BUF	1ST LOCATION OF ARRAY OF D.P. CHEBYSHEV COEFFICIENTS OF POSITION
C	    T   T(1) IS DP FRACTIONAL TIME IN INTERVAL COVERED BY
C		COEFFICIENTS AT WHICH INTERPOLATION IS WANTED
C		(0 .LE. T(1) .LE. 1).  T(2) IS DP LENGTH OF WHOLE
C		INTERVAL IN INPUT TIME UNITS.
C	NCF	# OF COEFFICIENTS PER COMPONENT
C	NCM	# OF COMPONENTS PER SET OF COEFFICIENTS (3 or 2)
C	NA	# OF SETS OF COEFFICIENTS IN FULL ARRAY
C		(I.E., # OF SUB-INTERVALS IN FULL INTERVAL)
C	IFL  INTEGER FLAG: =1 FOR POSITIONS ONLY
C			   =2 FOR POS AND VEL
C OUTPUTS:
C	pv	INTERPOLATED QUANTITIES REQUESTED.  DIMENSION
C		EXPECTED IS pv(NCM,IFL), DP.
C MODIFICATION HISTORY:
C-
	implicit double precision (A-H,O-Z)

	double precision BUF(NCF,NCM,*)
	double precision T(2)
	double precision pv(NCM,*)
	double precision PC(  18)	/1.D0,0.D0,16*0.D0/
	double precision VC(2:18)	/     1.D0,16*0.D0/
	double precision TWOT		/0.D0/
	integer		NP		/2/
	integer		NV		/3/
	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 = dble(NA)
	DT1 = dint(T(1))
	TC  = DNA*T(1)
	L   = idint(TC-DT1)+1

	TC = 2.D0*(dmod(TC,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(2) IS THE VALUE OF T1(TC) AND HENCE
	! CONTAINS THE VALUE OF TC ON THE PREVIOUS CALL.)

	if (TC .ne. PC(2)) then
	    TWOT = TC+TC
	    NP = 2
	    NV = 3
	    PC(2)= TC
	    VC(3)= TWOT+TWOT
	end if

	do I=NP+1,NCF			! At least NCF polynomials are needed
	    PC(I) = TWOT*PC(I-1)-PC(I-2)
	end do
	NP = max(NP,NCF)

	do I=1,NCM			! Interpolate positions
	    pv(I,1) = 0.D0
	    do J=NCF,1,-1
		pv(I,1) = pv(I,1)+PC(J)*BUF(J,I,L)
	    end do
	end do

	if (IFL .eq. 2) then
	    do I=NV+1,NCF		! At least NCF polynomials are needed
		VC(I) = TWOT*VC(I-1)+PC(I-1)+PC(I-1)-VC(I-2)
	    end do
	    NV = max(NV,NCF)

	    do I=1,NCM			! Interpolate velocities
		pv(I,2) = 0.D0
		do J=NCF,2,-1
		    pv(I,2) = pv(I,2)+VC(J)*BUF(J,I,L)
		end do
		pv(I,2) = pv(I,2)*((DNA+DNA)/T(2))
	    end do
	end if

	return
	end