/*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 "spquad.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	ONE	1.0e0
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
float /*FUNCTION*/ spquad(
long korder,
long npc,
float xi[],
float *pc,
float x1,
float x2)
{
#define PC(I_,J_)	(*(pc+(I_)*(korder)+(J_)))
	long int ii, im, left, mode;
	float a, aa, b, bb, denom, dx, q, s, spquad_v, ss[2], x;
	static long il1 = 1;
	static long il2 = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Ss = &ss[0] - 1;
	float *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 SPQUAD Krogh  Changes to use M77CON
	 *>> 1992-11-12 SPQUAD C. L. Lawson, JPL  Saving IL1 and IL2.
	 *>> 1992-10-27 C. L. Lawson, JPL
	 *>> 1988-03-16 C. L. Lawson, JPL
	 *
	 *   This subroutine computes the integral over [X1,X2] of a
	 *   piecewise polynomial using the piecewise Power
	 *   representation given by [XI, PC, NPC, KORDER].
	 *   The degrees of the polynomial pieces are KORDER-1.
	 *   Require  2 .le. KORDER .le. 20.
	 *   Permits X1 .le. X2 or X1 gt. X2.
	 *   The "proper interpolation interval" is from XI(1) to XI(NPC+1).
	 *   For most reliable results, X1 and X2 should lie in this
	 *   interval, however this subr will give a result even if this
	 *   is not the case by use of extrapolation.
	 *     ------------------------------------------------------------------
	 *  DESCRIPTION OF ARGUMENTS
	 *  INPUT:
	 *    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(,).
	 *           Require  2 .le. KORDER .le. 20.
	 *    NPC    - NUMBER OF POLYNOMIAL PIECES
	 *    XI()    - Breakpoint array of length NPC+1
	 *    PC(i,j) - Coeffs for the Power representation of a piecewise
	 *              polynomial.   i=1,KORDER , j=1,NPC
	 *              The coeffs PC(*,j) are used at XI(j) and on the
	 *              interval between XI(j) and XI(j+1).
	 *              PC(i,j) is the coefficient to be multiplied times
	 *              (X - X(j))**(i-1).
	 *    X1,X2  - END POINTS OF QUADRATURE INTERVAL, NORMALLY IN
	 *                XI(1) .LE. X .LE. XI(NPC+1) but extrapolation is
	 *                permitted.
	 *
	 *  OUTPUT:
	 *     SPQUAD - Integral of the piecewise polynomial from X1 to X2.
	 *     ------------------------------------------------------------------
	 *     Adapted from subroutine PPQUAD due to D. E. Amos,
	 *     Sandia, June, 1979.  Documented in SAND79-1825.  PPQUAD uses the
	 *     Taylor basis, whereas this subprogram uses the Power basis.
	 *     Uses the [XI, PC, NPC, KORDER] representation of a piecewise
	 *     polynomial as presented by Carl De Boor in
	 *     A PRACTICAL GUIDE TO SPLINES, Springer-Verlag, 1978.
	 *     (Called [XI, C, LXI, K] in the book.)
	 *     ------------------------------------------------------------------
	 *--S replaces "?": ?PQUAD, ?SFIND
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------
	 * */
	aa = fminf( x1, x2 );
	bb = fmaxf( x1, x2 );
	q = ZERO;
	if (aa == bb)
		goto L_90;
	ssfind( xi, 1, npc + 1, aa, &il1, &mode );
	ssfind( xi, 1, npc + 1, bb, &il2, &mode );
	for (left = il1; left <= il2; left++)
	{
		if (left == il1)
		{
			a = aa;
		}
		else
		{
			a = Xi[left];
		}
		if (left == il2)
		{
			b = bb;
		}
		else
		{
			b = Xi[left + 1];
		}
		x = a;
		for (ii = 1; ii <= 2; ii++)
		{
			Ss[ii] = ZERO;
			dx = x - Xi[left];
			if (dx != ZERO)
			{
				denom = (float)( korder );
				s = PC(left - 1,korder - 1)/denom;
				for (im = korder - 1; im >= 1; im--)
				{
					denom -= ONE;
					s = s*dx + PC(left - 1,im - 1)/denom;
				}
				Ss[ii] = s*dx;
			}
			x = b;
		}
		q += Ss[2] - Ss[1];
	}
	if (x1 > x2)
		q = -q;
L_90:
	;
	spquad_v = q;
	return( spquad_v );
#undef	PC
} /* end of function */