/*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 #include "fcrt.h" #include "ssbasi.h" #include /* 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 */