/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:43 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dsbasd.h" #include /* PARAMETER translations */ #define KMAX 20 /* end of PARAMETER translations */ void /*FUNCTION*/ dsbasd( long korder, long left, double t[], double x, long ideriv, double bderiv[]) { long int i, i1, id, istart, it1, it2, j, jp1, k1; double deltal[KMAX], deltar[KMAX], flk1, saved, term; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Bderiv = &bderiv[0] - 1; double *const Deltal = &deltal[0] - 1; double *const Deltar = &deltar[0] - 1; double *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 DSBASD Krogh Changes to use M77CON *>> 1993-01-12 DSBASD 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. * ------------------------------------------------------------------ *--D replaces "?": ?SBASD * Both versions use IERM1, IERV1 * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ id = ideriv; if (id < 0) { ierm1( "DSBASD", 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( "DSBASD", 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 = (double)( 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 */