/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:56 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "dpval.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	ONE	1.0e0
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
double /*FUNCTION*/ dpval(
long korder,
long npc,
double xi[],
double *pc,
double x,
long ideriv)
{
#define PC(I_,J_)	(*(pc+(I_)*(korder)+(J_)))
	long int i, mode;
	double dpval_v, dx, fac, fac1, fac2, val;
	static long lefti = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Xi = &xi[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1994-10-20 DPVAL Krogh  Changes to use M77CON
	 *>> 1992-10-27 DPVAL C. L. Lawson, JPL  Saving LEFTI
	 *>> 1988-03-16 C. L. Lawson, JPL
	 *
	 *     Computes value at X of the derivative of order IDERIV of a
	 *     piecewise polynomial given by coeffs relative to the Power basis.
	 *     The piecewise polynomial is specified by the arguments
	 *     XI, PC, NPC, KORDER.
	 *     The "proper interpolation interval" for X is from XI(1) to
	 *     XI(NPC+1), but extrapolation will be used outside this interval.
	 *     Based on subroutine PPVALU given on pp. 89-90 of
	 *     A PRACTICAL GUIDE TO SPLINES by Carl De Boor, Springer-Verlag,
	 *     1978, however PPVALU uses the Taylor basis and this subr uses
	 *     the Power basis.
	 *     The cases of IDERIV = 0, 1, or 2 are coded individually for
	 *     best efficiency.  The cases of 3 .le. IDERIV .lt. KORDER are
	 *     treated in a general way.
	 *     The result is zero when IDERIV .ge. KORDER.
	 *     ------------------------------------------------------------------
	 *     KORDER   [in]  Order of the polynomial pieces.  This is one
	 *           greater than the degree of the pieces.  KORDER is also the
	 *           leading dimension of the array PC(,).
	 *     NPC   [in]  Number of polynomial pieces.
	 *     (XI(j), j = 1, ..., NPC+1)  [in]  Breakpoints defining the
	 *           subintervals over which the polynomial pieces are defined.
	 *     ((PC(i,j), i = 1, ..., KORDER), j = 1, ..., NPC)  [in]
	 *           PC(*,j) contains the coefficients to be used to evaluate
	 *           the polynomial piece between XI(j) and XI(j+1).
	 *           PC(i,j) is the coefficient to be multiplied times
	 *           (X - X(j))**(i-1).
	 *     X  [in]  Abcissa at which function is to be evaluated.
	 *     IDERIV   [in]  Order of derivative to be evaluated.
	 *           Require IDERIV .ge. 0.  IDERIV = 0 means to evaluate the
	 *           function itself.  Derivatives of order .ge. KORDER will
	 *           be zero.
	 *     DPVAL [out]    The computed value or derivative value.
	 *     ------------------------------------------------------------------
	 *--D replaces "?": ?PVAL, ?SFIND
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	dsfind( xi, 1, npc + 1, x, &lefti, &mode );
	dx = x - Xi[lefti];
 
	if (ideriv == 0)
	{
		val = PC(lefti - 1,korder - 1);
		for (i = korder - 1; i >= 1; i--)
		{
			val = val*dx + PC(lefti - 1,i - 1);
		}
	}
	else if (ideriv == 1)
	{
		fac1 = (double)( korder - 1 );
		val = fac1*PC(lefti - 1,korder - 1);
		for (i = korder - 1; i >= 2; i--)
		{
			fac1 -= ONE;
			val = val*dx + fac1*PC(lefti - 1,i - 1);
		}
	}
	else if (ideriv == 2)
	{
		fac1 = (double)( korder - 1 );
		fac2 = fac1 - ONE;
		val = fac1*fac2*PC(lefti - 1,korder - 1);
		for (i = korder - 1; i >= 3; i--)
		{
			fac1 = fac2;
			fac2 -= ONE;
			val = val*dx + fac1*fac2*PC(lefti - 1,i - 1);
		}
	}
	else if (ideriv >= korder)
	{
		val = ZERO;
	}
	else
	{
		/*                           Here when 3 .le. IDERIV .lt. KORDER
		 * */
		fac1 = (double)( korder - 1 );
		fac2 = fac1;
		fac = fac1;
		for (i = 2; i <= ideriv; i++)
		{
			fac2 -= ONE;
			fac *= fac2;
		}
		val = fac*PC(lefti - 1,korder - 1);
		for (i = korder - 1; i >= (ideriv + 1); i--)
		{
			fac2 -= ONE;
			fac = (fac/fac1)*fac2;
			fac1 -= ONE;
			val = val*dx + fac*PC(lefti - 1,i - 1);
		}
	}
	dpval_v = val;
	return( dpval_v );
#undef	PC
} /* end of function */