;+
; NAME:
;	jpl_interp
; PURPOSE:
; 	Differentiates and interpolates a set of Chebyshev coefficients to give
;	position and velocity
; CATEGORY:
;	smei/gen/idl/ephem; JPL Lunar and Planetary Ephemeris
; CALLING SEQUENCE:
	FUNCTION jpl_interp, BUF,T,NCF,NCM,NA,IFL,PV
; INPUTS:
;	BUF 	1st location of array of double Chebyshev coefficients of position
;	T		T[1] is double fractional time in interval covered by
;			coefficients at which interpolation is wanted (0 <= T[1] <= 1).
;			T[2] is double length of whole interval in input time units
;	NCF		# coefficients per component
;	NCM		# components per set of coefficients
;	NA		# sets of coefficients in full array
;			(i.e, # sub-intervals in full interval)
;	IFL  	integer flag:	=1 for positions only
;							=2 for position and velocity
; OUTPUTS:
;	PV		array[NCM,IFL]; type: double
;				interpolated quantities requested.
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; COMMON BLOCKS:
	common	JPL_INTERP_SAVE, PC, VC, TWOT, NP, NV
; PROCEDURE:
; MODIFICATION HISTORY:
;	AUG-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		adapted from JPL Fortran software
;-

; 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