/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:48 */
/*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 "ssbasd.h"
#include <stdlib.h>
/* PARAMETER translations */
#define KMAX 20
/* end of PARAMETER translations */
void /*FUNCTION*/ ssbasd(
long korder,
long left,
float t[],
float x,
long ideriv,
float bderiv[])
{
long int i, i1, id, istart, it1, it2, j, jp1, k1;
float deltal[KMAX], deltar[KMAX], flk1, saved, term;
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const Bderiv = &bderiv[0] - 1;
float *const Deltal = &deltal[0] - 1;
float *const Deltar = &deltar[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 SSBASD Krogh Changes to use M77CON
*>> 1993-01-12 SSBASD CLL Bringing all computation inline.
*>> 1992-11-02 C. L. Lawson, JPL
*>> 1988-03-22 C. L. Lawson, JPL
*
* Given the knot-set T(i), i = 1, ..., NT, this subr computes
* values at X of the derivative of order IDERIV of each of the
* the KORDER B-spline basis functions of order KORDER
* that are nonzero on the open interval (T(LEFT), T(LEFT+1).
* Require T(LEFT) < T(LEFT+1) and KORDER .le. LEFT .le.
* NT-KORDER.
* For proper evaluation, X must lie in the closed interval
* [T(LEFT), T(LEFT+1], otherwise the computation constitutes
* extrapolation.
* The basis functions whose derivatives are returned will be those
* indexed from LEFT+1-KORDER through LEFT, and their values will be
* stored in BDERIV(I), I = 1, ..., KORDER.
* ------------------------------------------------------------------
* Method:
*
* In general there are two stages: First compute values of basis
* functions of order KORDER-IDERIV. Then, if IDERIV > 0, transform
* these values to produce values of higher order derivatives of
* higher order basis functions, ultimately obtaining values of
* the IDERIV derivative of basis functions of order KORDER.
*
* The first stage uses the recursion:
*
* X - T(I) T(I+J+1) - X
* B(I,J+1)(X) = -----------B(I,J)(X) + ---------------B(I+1,J)(X)
* T(I+J)-T(I) T(I+J+1)-T(I+1)
*
* where B(I,J) denotes the basis function of order J and index I.
* For order J, the only basis functions that can be nonzero on
* [T(LEFT), T(LEFT+1] are those indexed from LEFT+1-J to LEFT.
* For order 1, the only basis function nonzero on this interval is
* is B(LEFT,1), whose value is 1. The organization of the calculation
* using the above recursion follows Algorithm (8) in Chapter X of the
* book by DeBoor.
*
* For the second stage, let B(ID, K, J) denote the value at X of the
* IDth derivative of the basis function of order K, with index J.
* From the first stage we will have values of B(0, KORDER-IDERIV, j)
* for j = LEFT+1-KORDER+IDERIV), ..., LEFT, stored in BDERIV(i), for
* i = 1+IDERIV, ..., KORDER.
*
* Loop for ID = 1, ..., IDERIV
* Replace the contents of BDERIV() by values of
* B(ID, KORDER-IDERIV+ID, j) for
* j = LEFT+1-KORDER+IDERIV-ID), ..., LEFT, storing these in
* BDERIV(i), i = 1+IDERIV-ID, ..., KORDER.
* End loop
*
* The above loop uses formula (10), p. 138, of
* A PRACTICAL GUIDE TO SPLINES by Carl DeBoor, Springer-Verlag,
* 1978, when ID = 1, and successive derivatives of that formula
* when ID > 1. Note that we are using Formula (10) in the special
* case in which (Alpha sub i) = 1, and all other Alpha's are zero.
*
* This approach and implementation by C. L. Lawson, JPL, March 1988.
* ------------------------------------------------------------------
* KORDER [in] Gives both the order of the B-spline basis functions
* whose derivatives are to be evaluated,
* and the number of such functions to be evaluated.
* Note that the polynomial degree of these functions is one less
* than the order. Example: For cubic splines the order is 4.
* Require 1 .le. KORDER .le. KMAX, where KMAX is an internal
* parameter.
* LEFT [in] Index identifying the left end of the interval
* [T(LEFT), T(LEFT+1)] relative to which the B-spline basis
* functions will be evaluated.
* Require KORDER .le. LEFT .le. NT-KORDER
* and T(LEFT) < T(LEFT+1)
* The evaluation is proper if T(LEFT) .le. X .le. T(LEFT+1), and
* otherwise constitutes extrapolation.
* DIVISION BY ZERO will result if T(LEFT) = T(LEFT+1)
* T() [in] Knot sequence, indexed from 1 to NT, with
* NT .ge. LEFT+KORDER. Knot values must be nonincreasing.
* Repetition of values is permitted and has the effect of
* decreasing the order of contimuity at the repeated knot.
* Proper function representation by splines of order K is
* supported on the interval from T(K) to T(NT+1-K).
* Extrapolation can be done outside this interval.
* X [in] The abcissa at which the B-spline basis functions are to
* be evaluated. The evaluation is proper if T(KORDER) .le. X .le
* T(NT+1-KORDER), and constitutes extrapolation otherwise.
* IDERIV [in] Order of derivative requested. Require IDERIV .ge. 0.
* Derivatives of order .ge. KORDER are zero.
* BDERIV() [out] On normal return, this array will contain in
* BDERIV(i), i = 1, ...,KORDER, the values computed by this subr.
* These are values at X of the derivative of order IDERIV of each
* of the basis functions indexed LEFT+1-KORDER through LEFT.
* ------------------------------------------------------------------
*--S replaces "?": ?SBASD
* Both versions use IERM1, IERV1
* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
id = ideriv;
if (id < 0)
{
ierm1( "SSBASD", 2, 2, "Require IDERIV .ge. 0", "IDERIV",
id, '.' );
return;
}
if (id >= korder)
{
for (i = 1; i <= korder; i++)
{
Bderiv[i] = 0.0e0;
}
return;
}
if (korder > KMAX)
{
ierm1( "SSBASD", 1, 2, "Require KORDER .le. KMAX.", "KORDER"
, korder, ',' );
ierv1( "KMAX", KMAX, '.' );
return;
}
istart = 1 + id;
k1 = korder - id;
/* Evaluate K1 basis functions of order K1. Store values in
* BDERIV(ISTART:ISTART+K1-1)
* */
Bderiv[istart] = 1.0e0;
for (j = 1; j <= (k1 - 1); j++)
{
jp1 = j + 1;
Deltar[j] = T[left + j] - x;
Deltal[j] = x - T[left + 1 - j];
saved = 0.0e0;
for (i = 1; i <= j; i++)
{
term = Bderiv[id + i]/(Deltar[i] + Deltal[jp1 - i]);
Bderiv[id + i] = saved + Deltar[i]*term;
saved = Deltal[jp1 - i]*term;
}
Bderiv[id + jp1] = saved;
}
/* Loop IDERIV times, each time advancing the order of the
* derivative by one, so final results are values of derivatives
* of order IDERIV.
* */
for (i1 = istart; i1 >= 2; i1--)
{
flk1 = (float)( k1 );
it1 = left + 1 - i1;
it2 = it1 - k1;
for (i = i1; i <= korder; i++)
{
Bderiv[i] = Bderiv[i]*flk1/(T[it1 + i] - T[it2 + i]);
}
Bderiv[i1 - 1] = -Bderiv[i1];
for (i = i1; i <= (korder - 1); i++)
{
Bderiv[i] -= Bderiv[i + 1];
}
k1 += 1;
}
return;
} /* end of function */