#  INTRP_STATE -- state vector interpolation for a specified epoch
#
#  Description:
#  ------------
#  Using Lagrange's formula of polynomial interpolation, obtain the state
#  vector of any epoch from the input state vector table.
#  The input state vectors must be tabulated with equal epoch interval.
#
#  Date 	Author 		Description
#  ----		------		-----------
#  28-Feb-1988  J.-C. Hsu       original, vxyzin.for
#  17-Apr-1990	J.-C. Hsu	rewrite in SPP   
#-------------------------------------------------------------------------------

procedure intrp_state (epoch, x, y, z, first, grid, xyz, npts, nlag)

double	epoch		# input: desired epoch
double	x[npts], y[npts], z[npts]
			# input: state vectors
double	first		# input: epoch of the first state vector
double	grid		# input: epoch increment of the state vector array
double	xyz[3]		# output: state vector of the desired epoch
int	npts		# input: number of entries of the input state vectors
int	nlag		# input: number of points used in Lagrange's
			#        polynomial interpolation formula

int	j, icenter, istart, istop
double	dxyz[3]		# error of xyz
pointer	sp, ta, c, d
#------------------------------------------------------------------------------
begin

	# allocate space for the working arrays
	call smark (sp)
	call salloc (ta, nlag, TY_DOUBLE)
	call salloc (c, nlag, TY_DOUBLE)
	call salloc (d, nlag, TY_DOUBLE)

	# find the array index of the desired epoch
	icenter = int ((epoch - first) / grid)
        istart = icenter - nlag / 2
        istop = istart + nlag - 1

	# make sure the starting and stopping points are within limits 
        if (istart < 1 || istop > npts)
	    call error (1, "input epoch out of ephemeris range")

	# construct the epoch array used in interpolation
        do j = 0, nlag-1
            Memd[ta+j] = first + grid * (istart-1+j)

	# perform the interpolation
        call vpolin (Memd[ta], x[istart], nlag, Memd[c], Memd[d], epoch,    
			xyz[1], dxyz[1])
        call vpolin (Memd[ta], y[istart], nlag, Memd[c], Memd[d], epoch,    
			xyz[2], dxyz[2])
        call vpolin (Memd[ta], z[istart], nlag, Memd[c], Memd[d], epoch,    
			xyz[3], dxyz[3])

	call sfree (sp)
end