/*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 "ssbasi.h"
#include <stdlib.h>
/* PARAMETER translations */
#define HALF 0.50e0
#define KMAX 20
#define ZERO 0.0e0
/* end of PARAMETER translations */
void /*FUNCTION*/ ssbasi(
long k,
long n,
float t[],
float x1,
float x2,
long *j1,
long *j2,
float basi[])
{
long int ik, il1, il2, j, jf, left, m, mf, mode;
float a, aa, b, bb, bma, bpa, bval1[KMAX], bval2[KMAX], c1, fac,
ta, tb;
static float gpts[9]={0.5773502691896257645091487805019574556476e0,
0.2386191860831969086305017216807119354186e0,0.6612093864662645136613995950199053470064e0,
0.9324695142031520278123015544939946091348e0,0.1488743389816312108848260011297199846176e0,
0.4333953941292471907992659431657841622001e0,0.6794095682990244062343273651148735757693e0,
0.8650633666889845107320966884234930485275e0,0.9739065285171717200779640120844520534283e0};
static float gwts[9]={1.0e0,0.4679139345726910473898703439895509948120e0,
0.3607615730481386075698335138377161116612e0,0.1713244923791703450402961421727328935273e0,
0.2955242247147528701738929946513383294210e0,0.2692667193099963550912269215694693528586e0,
0.2190863625159820439955349342281631924634e0,0.1494513491505805931457763396576973323908e0,
0.06667134430868813759356880989333179285686e0};
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const Basi = &basi[0] - 1;
float *const Bval1 = &bval1[0] - 1;
float *const Bval2 = &bval2[0] - 1;
float *const Gpts = &gpts[0] - 1;
float *const Gwts = &gwts[0] - 1;
float *const T = &t[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 SSBASI Krogh Changes to use M77CON
*>> 1993-01-12 SSBASI CLL Using SSBASD in place of _SBAS.
*>> 1992-11-12 C. L. Lawson, JPL Made J1 & J2 inout. Were out.
*>> 1992-11-02 C. L. Lawson, JPL
*>> 1988-04-01 C. L. Lawson, JPL
*
* This subroutine computes the integral over [X1,X2] of each of the N
* members of a Kth order family of B-spline basis functions defined
* over the knot set T(). The values of these integrals will be
* returned in BASI(j), j = 1, ..., N. The indices J1 and J2 will be
* set to indicate the subset of these that might be nonzero. Thus
* BASI(j), j = J1, ..., J2 might be nonzero, whereas the other
* elements of BASI() will be zero.
* Orders K as high as 20 are permitted.
* Permits X1 .le. X2 or X1 gt. X2. In the former case all returned
* values will be .ge. 0, while in the latter case all will be .le. 0.
* The "proper interpolation interval" is from T(K) to T(N+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.
* This subr applies a 2, 6, or 10 point Gauss formula on
* subintervals of [X1,X2] which are formed by the included
* (distinct) knots.
* ------------------------------------------------------------------
* ARGUMENTS
*
* K [in] ORDER OF B-SPLINE, 2 .le. K .le. 20
* Note that the polynomial degree of the segments of
* the spline are one less than the order.
* N [in] LENGTH OF COEFFICIENT ARRAY
* T() [in] KNOT ARRAY OF LENGTH N+K
* X1,X2 [in] END POINTS OF QUADRATURE INTERVAL.
* Permit X1 .le. X2 or X1 .gt. X2.
*
* J1, J2 [inout integers] On entry must contain integer values.
* Will be passed to SSFIND for lookup. Any integer values
* are ok but previous returned values may aid efficiency.
* On return will satisfy 1 .le. J1 .le. J2 .le. N.
* Indicates that only the output values BASI(j), j = J1,...,J2
* might be nonzero.
* BASI() [out, floating pt] Array in which the N integrals will be
* stored.
*
* ERROR CONDITIONS: None
* ------------------------------------------------------------------
* Based on spline quadrature subroutines called BQUAD or BSQUAD
* written by D. E. Amos, Sandia, June, 1979. Documented in
* SAND79-1825.
* Uses the (T, BCOEF, N, K) representation of a B-spline as
* presented in A PRACTICAL GUIDE TO SPLINES by Carl De Boor,
* Springer-Verlag, 1978.
* Current version by C. L. Lawson, JPL, March 1988.
* ------------------------------------------------------------------
*--S replaces "?": ?SBASI, ?SFIND, ?SBASD
* ------------------------------------------------------------------
* GPTS(), GWTS() Index 1 selects 2-point Gaussian formula.
* Indices 2-4 select 6-point Gaussian formula.
* Indices 5-9 select 10-point Gaussian formula.
* ------------------------------------------------------------------ */
/* DATA GPTS / 5.77350269189626E-01, 2.38619186083197E-01,
* 1 6.61209386466265E-01, 9.32469514203152E-01, 1.48874338981631E-01,
* 2 4.33395394129247E-01, 6.79409568299024E-01, 8.65063366688985E-01,
* 3 9.73906528517172E-01/
*
* DATA GWTS / 1.00000000000000E+00, 4.67913934572691E-01,
* 1 3.60761573048139E-01, 1.71324492379170E-01, 2.95524224714753E-01,
* 2 2.69266719309996E-01, 2.19086362515982E-01, 1.49451349150581E-01,
* 3 6.66713443086880E-02/
*
* Gaussian quadrature abcissae and weights to 40 decimal digits. */
/* ------------------------------------------------------------------ */
if (k > KMAX)
{
ierm1( "SSBASI", 1, 2, "Require KORDER .le. KMAX.", "KORDER"
, k, ',' );
ierv1( "KMAX", KMAX, '.' );
return;
}
aa = fminf( x1, x2 );
bb = fmaxf( x1, x2 );
il1 = *j1 + k - 1;
il2 = *j2;
ssfind( t, k, n + 1, aa, &il1, &mode );
ssfind( t, k, n + 1, bb, &il2, &mode );
if ((bb == T[il2] && il2 > k) && il2 > il1)
il2 -= 1;
*j1 = il1 + 1 - k;
*j2 = il2;
for (j = 1; j <= n; j++)
{
Basi[j] = ZERO;
}
if (aa == bb)
return;
/* Selection of Gaussian quadrature formula:
* KORDER .le. 4 => 2 point formula
* 5 .le. KORDER .le. 12 => 6 point formula
* KORDER .ge. 13 => 10 point formula
* */
if (k <= 4)
{
jf = 0;
mf = 1;
}
else if (k <= 12)
{
jf = 1;
mf = 3;
}
else
{
jf = 4;
mf = 5;
}
for (left = il1; left <= il2; left++)
{
ta = T[left];
tb = T[left + 1];
if (ta != tb)
{
if (left == k)
{
a = aa;
}
else
{
a = fmaxf( aa, ta );
}
if (left == n)
{
b = bb;
}
else
{
b = fminf( bb, tb );
}
bma = HALF*(b - a);
bpa = HALF*(b + a);
for (m = 1; m <= mf; m++)
{
c1 = bma*Gpts[jf + m];
fac = Gwts[jf + m]*bma;
ssbasd( k, left, t, bpa + c1, 0, bval1 );
ssbasd( k, left, t, bpa - c1, 0, bval2 );
j = left - k;
for (ik = 1; ik <= k; ik++)
{
j += 1;
Basi[j] += fac*(Bval1[ik] + Bval2[ik]);
}
}
}
}
if (x1 > x2)
{
for (j = *j1; j <= *j2; j++)
{
Basi[j] = -Basi[j];
}
}
return;
} /* end of function */