/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:59 */
/*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 "ssvala.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	KMAX	20
#define	ONE	1.0e0
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ ssvala(
long k,
long n,
float t[],
long nderiv,
float *bdif,
float x,
float svalue[])
{
#define BDIF(I_,J_)	(*(bdif+(I_)*(n)+(J_)))
	long int i, id, j, jp1, kp1mn, l, left, leftpl, mode, nder;
	float deltal[KMAX], deltar[KMAX], saved, term, vnikx[KMAX];
	static long lefti = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Deltal = &deltal[0] - 1;
	float *const Deltar = &deltar[0] - 1;
	float *const Svalue = &svalue[0] - 1;
	float *const T = &t[0] - 1;
	float *const Vnikx = &vnikx[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 SSVALA Krogh  Changes to use M77CON
	 *>> 1993-01-12 SSVALA CLL Bringing basis computation inline.
	 *>> 1992-11-12 C. L. Lawson, JPL Saving LEFTI.
	 *>> 1992-11-02 C. L. Lawson, JPL
	 *>> 1988-03-16 C. L. Lawson, JPL
	 *
	 *     A spline function is defined by the contents of K, N, and T().
	 *     Given array BDIF(,) containing B-spline coefficients and
	 *     differences of these, this subroutine computes and stores into
	 *     SVALUE(1:NDERIV+1) the value and derivatives through order NDERIV
	 *     of the spline function evaluated at the abcissa, X.
	 *     Derivatives of order .ge. K will be zero.
	 *
	 *     The subroutine, BSPLPP, given on pp. 140-141 of
	 *     A PRACTICAL GUIDE TO SPLINES by Carl De Boor, Springer-Verlag,
	 *     1978, has been recoded as separate subroutines: DSTOT (or DSTOP)
	 *     calling DSDIF and SSVALA.  This subroutine has the functionality
	 *     of lines 72-103 of BSPLPP in the book.
	 *     ------------------------------------------------------------------
	 *--S replaces "?": ?SVALA, ?SFIND
	 *     Both versions use IERM1, IERV1
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	if (k > KMAX)
	{
		ierm1( "SSVALA", 1, 2, "Require KORDER .le. KMAX.", "KORDER"
		 , k, ',' );
		ierv1( "KMAX", KMAX, '.' );
		return;
	}
	for (i = 1; i <= (nderiv + 1); i++)
	{
		Svalue[i] = ZERO;
	}
	ssfind( t, k, n + 1, x, &lefti, &mode );
 
	/*     LEFTI has been found so that [T(LEFTI), T(LEFTI+1)] is the
	 *     reference interval for the polynomial piece to be used in
	 *     evaluating at X.
	 * */
	nder = min( nderiv, k - 1 );
	kp1mn = k - nder;
 
	/*        Compute values of KP1MN basis functions of order KP1MN at X.
	 *        Store these in VNIKX(1:KP1MN).  Here KP1MN .ge. 1.
	 * */
	Vnikx[1] = ONE;
	for (j = 1; j <= (kp1mn - 1); j++)
	{
		jp1 = j + 1;
		Deltar[j] = T[lefti + j] - x;
		Deltal[j] = x - T[lefti + 1 - j];
		saved = ZERO;
		for (i = 1; i <= j; i++)
		{
			term = Vnikx[i]/(Deltar[i] + Deltal[jp1 - i]);
			Vnikx[i] = saved + Deltar[i]*term;
			saved = Deltal[jp1 - i]*term;
		}
		Vnikx[jp1] = saved;
	}
 
	/*        Loop on ID.  For each value of ID, evaluate the derivative of
	 *        order ID-1 at X, and store it in SVALUE(ID).
	 *        Then, if ID > 1, update basis function values
	 *        in VNIKX() to the next higher order. (From order J to JP1).
	 * */
	for (id = nder + 1; id >= 1; id--)
	{
		left = lefti - kp1mn;
		for (l = 1; l <= kp1mn; l++)
		{
			leftpl = left + l;
			Svalue[id] += Vnikx[l]*BDIF(id - 1,leftpl - 1);
		}
		if (id > 1)
		{
			j = kp1mn;
			jp1 = j + 1;
			Deltar[j] = T[lefti + j] - x;
			Deltal[j] = x - T[lefti + 1 - j];
			saved = ZERO;
			for (i = 1; i <= j; i++)
			{
				term = Vnikx[i]/(Deltar[i] + Deltal[jp1 - i]);
				Vnikx[i] = saved + Deltar[i]*term;
				saved = Deltal[jp1 - i]*term;
			}
			Vnikx[jp1] = saved;
		}
		kp1mn += 1;
	}
	return;
#undef	BDIF
} /* end of function */