/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:04 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dc2bas.h" #include /* PARAMETER translations */ #define C18 18.0e0 #define C36 36.0e0 #define C6 6.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dc2bas( double x, long nderiv, long iseg, LOGICAL32 *newseg, double b[], long nb, double p[]) { long int i, i2, ib, ib1, ib2, is, ism4, iu, j, mb, ndp1; double cs, flnseg, h, xs; static double scale, u[8]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const B = &b[0] - 1; double *const P = &p[0] - 1; double *const U = &u[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-11-11 DC2BAS Krogh Declared all vars. *>> 1994-10-20 DC2BAS Krogh Changes to use M77CON *>> 1987-12-09 DC2BAS Lawson Initial code. *--D replaces "?": ?C2BAS * * Given X in the segment ISEG, this subr evaluates at X the * four cubic B-spline basis functions that are not identically * zero throughout segment ISEG. * This is a reorganization of a set of three subroutines * SBAS1, SBAS2, & SBAS3 used in LIB*JPL$ from about 1970 on. * The reorganization improves portability and brings * conformance to Fortran 77. * More features are avalable from a similar subr in a package * due to Carl de Boor, however this subr is being preserved * for convenience of usage with some other JPL subrs -- * particularly the size and interpretation of the B() array. * C. L. Lawson, JPL, 1984 July 9, July 1987. * ------------------------------------------------------------------ * A R G U M E N T S * * X [in] Argument at which basis fcns are to be evaluated. * Must satisfy B(ISEG) .le. X .le. B(ISEG+1) * * NDERIV [in] = 0, 1, 2, or 3. Indicates the order of derivative * to be computed. 0 selects function value, * 1 selects 1st derivative, etc. * * ISEG [in] Index of current segment. Require 1 .le. ISEG .le. NB-1 * Segment ISEG is the interval [B(ISEG), B(ISEG+1)]. * The current X must be in this segment. * * NEWSEG [inout] = .true. means ISEG is not the same as on * previous call. In this case this subr will set new values * into U() and reset NEWSEG = .false. The user prog must * set NEWSEG = .true. on the first call to this subr, and * on any subsequent call at which the user prog has set * ISEG to a new value. * * = .false. means this subr can assume the values in U() * are correct for the current ISEG. * * B() [in] Array of knots, including the left and right * endpoints. Contains B(I), I = 1, NB. Set by user. * Not changed by this subr. Must be strictly increasing, * i.e. require B(i) .lt. B(i+1), i = 1, NB-1. * * NB [in] Number of knots, including endpoints. The number * of degrees of freedom of the set of cubic splines * defined over this set of knots is NB+2. * * P() [out] Array of length 4 into which this subr will store the * values of the 4 basis fcns that are not identically * zero over segment ISEG. The indices of these basis * fcns are ISEG, ISEG+1, ISEG+2, & ISEG+3. The values * will be stored in (P(i), i=1,4). * ------------------------------------------------------------------ * IMPORTANT INTERNAL VARIABLES * * U() Array of length 8 in which this subr saves a scaled * record of the knots near the current segment. * * SCALE Scale factor for abcissas. Computed as the reciprocal * of the average segment length. * ------------------------------------------------------------------ * */ /* ------------------------------------------------------------------ */ mb = nb; /* The following statement does a type conversion. */ flnseg = mb - 1; if (*newseg) { *newseg = FALSE; scale = flnseg/(B[nb] - B[1]); is = iseg; ism4 = is - 4; ib1 = max( 1, is - 3 ); ib2 = min( mb, is + 4 ); /* COPY PART OF B() INTO U() */ for (ib = ib1; ib <= ib2; ib++) { iu = ib - ism4; U[iu] = scale*B[ib]; } /* ADJOIN ARTIFICIAL SEGMENTS AT LEFT END */ iu = ib1 - ism4 - 1; if (iu != 0) { h = U[iu + 2] - U[iu + 1]; L_30: U[iu] = U[iu + 1] - h; iu -= 1; if (iu > 0) goto L_30; } /* ADJOIN ARTIFICIAL SEGMENTS AT RIGHT END */ iu = ib2 - ism4 + 1; if (iu <= 8) { h = U[iu - 1] - U[iu - 2]; L_50: U[iu] = U[iu - 1] + h; iu += 1; if (iu <= 8) goto L_50; } } /* ------------------------------------------------------------------ */ ndp1 = nderiv + 1; xs = x*scale; switch (ndp1) { case 1: goto L_70; case 2: goto L_90; case 3: goto L_110; case 4: goto L_130; } L_70: for (i = 1; i <= 4; i++) { P[i] = powi(xs - U[i],3); } cs = C6; goto L_150; L_90: for (i = 1; i <= 4; i++) { P[i] = powi(xs - U[i],2); } cs = C18*scale; goto L_150; L_110: for (i = 1; i <= 4; i++) { P[i] = xs - U[i]; } cs = C36*SQ(scale); goto L_150; L_130: for (i = 1; i <= 4; i++) { P[i] = C36; } cs = CUBE(scale); L_150: ; /* COMPUTE FOURTH DIVIDED DIFFERENCES * */ for (j = 1; j <= 4; j++) { for (i = 1; i <= 3; i++) { i2 = i + j; P[i] = (P[i + 1] - P[i])/(U[i2] - U[i]); } P[4] = -P[4]/(U[j + 4] - U[4]); } for (i = 1; i <= 4; i++) { P[i] *= cs; } return; } /* end of function */