/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:09 */
/*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 "shint.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	ONE	1.0e0
#define	SIX	6.0e0
#define	THREE	3.0e0
#define	TWO	2.0e0
		/* end of PARAMETER translations */
 
float /*FUNCTION*/ shint(
float x,
long nderiv,
long ntab,
float xtab[],
float ytab[],
float yptab[])
{
	long int look1;
	float q0, q1, q2, q3, q4, q5, r, shint_v;
	static long look = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Xtab = &xtab[0] - 1;
	float *const Yptab = &yptab[0] - 1;
	float *const Ytab = &ytab[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-11-11 SHINT  Krogh   Declared all vars.
	 *>> 1994-10-20 SHINT  Krogh  Changes to use M77CON
	 *>> 1987-12-09 SHINT  Lawson  Initial code.
	 *--S replaces "?": ?HINT
	 *
	 *          This subprogram does look-up and Hermite cubic interpolation
	 *     to evaluate a function at X defined by the data
	 *     NTAB, XTAB(), YTAB(), & YPTAB().
	 *     The values in XTAB() must be either strictly increasing or
	 *     strictly decreasing.
	 *     Optionally this subprogram will evaluate either the function or
	 *     its first, second, or third derivative.
	 *
	 *     The look-up method is a linear search beginning with the index,
	 *     LOOK, saved internally from the previous reference to this
	 *     subprogram.  LOOK is reset if it exceeds the current NTAB.
	 *     Evaluation at an interior tabular point uses data from the
	 *     interval in the direction of decreasing abcissa values.
	 *     Evaluation outside the table (extrapolation) uses data from the
	 *     nearest tabular interval.
	 *     ------------------------------------------------------------------
	 *     Based on subprograms FXDP and DYDXDP by Lawson, Hanson, Lang, and
	 *     Campbell, JPL, 1968-1974.
	 *     Modified July 1984 and July 1987 by C. L. Lawson, JPL, for
	 *     inclusion in the MATH 77 library.
	 *     ------------------------------------------------------------------
	 *     X  [in]  Argument at which interpolation is to be done.
	 *           For best results X should be between XTAB(1) and XTAB(NTAB),
	 *           however this subprogram will give an extrapolated
	 *           value when X is outside this range.
	 *
	 *     NDERIV  [in]   = 0 to compute function value.
	 *                    = 1 to compute value of first derivative.
	 *                    = 2 to compute value of second derivative.
	 *                    = 3 to compute value of third derivative.
	 *
	 *     NTAB  [in]  No. of points in XTAB(), YTAB(), and YPTAB().
	 *                 Require NTAB .ge. 2.
	 *
	 *     XTAB()  [in]  Array of abcissas.  Values must either be strictly
	 *                   increasing or strictly decreasing.
	 *
	 *     YTAB()  [in]  Array of function values.
	 *                   YTAB(i) = f(XTAB(i))
	 *
	 *     YPTAB()  [in]  Array of first derivative values.
	 *                    YPTAB(i) = f'(XTAB(i))
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------
	 *
	 * */
	if (Xtab[1] < Xtab[2])
	{
		/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		 *                             Here when values in XTAB() are increasing. */
		if (look >= ntab)
			look = max( 1, ntab/2 );
		/*                                      Here  1 .le. LOOK .le. NTAB-1 */
		if (x > Xtab[look + 1])
			goto L_30;
		/*                 Search toward decreasing abcissae, decreasing indices. */
L_25:
		if (look == 1)
			goto L_40;
		if (x > Xtab[look])
			goto L_40;
		look -= 1;
		goto L_25;
		/*                 Search toward increasing abcissae, increasing indices. */
L_30:
		if (look == ntab - 1)
			goto L_40;
L_35:
		look += 1;
		if (look == ntab - 1)
			goto L_40;
		if (x > Xtab[look + 1])
			goto L_35;
L_40:
		look1 = look + 1;
	}
	else
	{
		/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		 *                             Here when values in XTAB() are decreasing. */
		if (look > ntab || look < 2)
			look = max( 2, ntab/2 );
 
		/*                         Here  2 .le. LOOK .le. NTAB
		 * */
		if (x > Xtab[look - 1])
			goto L_50;
		/*                 Search toward decreasing abcissae, increasing indices. */
L_45:
		if (look == ntab)
			goto L_60;
		if (x > Xtab[look])
			goto L_60;
		look += 1;
		goto L_45;
		/*                 Search toward increasing abcissae, decreasing indices. */
L_50:
		if (look == 2)
			goto L_60;
L_55:
		look -= 1;
		if (look == 2)
			goto L_60;
		if (x > Xtab[look - 1])
			goto L_55;
L_60:
		look1 = look - 1;
	}
	/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	 *              Here XTAB(LOOK) < XTAB(LOOK1) and X is in this interval
	 *              unless we are extrapolating.
	 *              Evaluate the Hermite interpolation formulas.
	 * */
	r = Xtab[look1] - Xtab[look];
	q0 = (x - Xtab[look])/r;
	q1 = q0 - ONE;
 
	switch (nderiv + 1)
	{
		case 1: goto L_100;
		case 2: goto L_101;
		case 3: goto L_102;
		case 4: goto L_103;
	}
 
L_100:
	q4 = TWO*q0;
	q2 = ONE + q4;
	q3 = THREE - q4;
	q4 = SQ(q1);
	q5 = SQ(q0);
	shint_v = q4*(q2*Ytab[look] + r*q0*Yptab[look]) + q5*(q3*Ytab[look1] +
	 r*q1*Yptab[look1]);
	return( shint_v );
 
L_101:
	q2 = THREE*q0;
	q3 = SIX*q0/r;
	shint_v = q1*(q3*(Ytab[look] - Ytab[look1]) + (q2 - ONE)*Yptab[look]) +
	 q0*(q2 - TWO)*Yptab[look1];
	return( shint_v );
 
L_102:
	q2 = q1 + q0;
	shint_v = (TWO/r)*((THREE/r)*(Ytab[look] - Ytab[look1])*q2 + Yptab[look]*
	 (q2 + q1) + Yptab[look1]*(q2 + q0));
 
	return( shint_v );
 
L_103:
	shint_v = (SIX/SQ(r))*((TWO/r)*(Ytab[look] - Ytab[look1]) + Yptab[look] +
	 Yptab[look1]);
	return( shint_v );
 
} /* end of function */