/* lgresp.f -- translated by f2c (version 19980913). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* $Procedure LGRESP ( Lagrange interpolation on equally spaced points ) */ doublereal lgresp_(integer *n, doublereal *first, doublereal *step, doublereal *yvals, doublereal *work, doublereal *x) { /* System generated locals */ integer i__1, i__2; doublereal ret_val; /* Local variables */ doublereal newx; integer i__, j; extern /* Subroutine */ int chkin_(char *, ftnlen); doublereal c1, c2; extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *, ftnlen), setmsg_(char *, ftnlen), errint_(char *, integer *, ftnlen); extern logical return_(void); /* $ Abstract */ /* Evaluate a Lagrange interpolating polynomial for a specified */ /* set of coordinate pairs whose first components are equally */ /* spaced, at a specified abcissisa value. */ /* $ Copyright */ /* Copyright (1995), California Institute of Technology. */ /* U.S. Government sponsorship acknowledged. */ /* $ Required_Reading */ /* None. */ /* $ Keywords */ /* INTERPOLATION */ /* POLYNOMIAL */ /* $ Declarations */ /* $ Brief_I/O */ /* Variable I/O Description */ /* -------- --- -------------------------------------------------- */ /* N I Number of points defining the polynomial. */ /* FIRST I First abscissa value. */ /* STEP I Step size. */ /* YVALS I Ordinate values. */ /* WORK I-O Work space array. */ /* X I Point at which to interpolate the polynomial. */ /* The function returns the value at X of the unique polynomial of */ /* degree N-1 that fits the points in the plane defined by FIRST, */ /* STEP, and YVALS. */ /* $ Detailed_Input */ /* N is the number of points defining the polynomial. */ /* The arrays XVALS and YVALS contain N elements. */ /* FIRST, */ /* STEP are, respectively, a starting abscissa value and a */ /* step size that define the set of abscissa values */ /* at which a Lagrange interpolating polynomial is to */ /* be defined. The set of abscissa values is */ /* FIRST + I * STEP, I = 0, ..., N-1 */ /* STEP must be non-zero. */ /* YVALS is an array of ordinate values that, together with */ /* the abscissa values defined by FIRST and STEP, */ /* define N ordered pairs belonging to the graph of */ /* a function. The set of points */ /* ( FIRST + (I-1)*STEP, YVALS(I) ) */ /* where I ranges from 1 to N, define the Lagrange */ /* polynomial used for interpolation. */ /* WORK is a work space array of the same dimension as */ /* XVALS and YVALS. It is used by this routine as a */ /* scratch area to hold intermediate results. WORK */ /* is permitted to coincide with YVALS. */ /* X is the abscissa value at which the interpolating */ /* polynomial is to be evaluated. */ /* $ Detailed_Output */ /* The function returns the value at X of the unique polynomial of */ /* degree N-1 that fits the points in the plane defined by FIRST, */ /* STEP, and YVALS. */ /* $ Parameters */ /* None. */ /* $ Exceptions */ /* 1) If STEP is zero, the error SPICE(INVALIDSTEPSIZE) will */ /* be signalled. The function will return the value 0.D0. */ /* 2) If N is less than 1, the error SPICE(INVALIDSIZE) is */ /* signalled. The function will return the value 0.D0. */ /* 3) This routine does not attempt to ward off or diagnose */ /* arithmetic overflows. */ /* $ Files */ /* None. */ /* $ Particulars */ /* Given a set of N distinct abscissa values and corresponding */ /* ordinate values, there is a unique polynomial of degree N-1, */ /* often called the `Lagrange polynomial', that fits the graph */ /* defined by these values. The Lagrange polynomial can be used to */ /* interpolate the value of a function at a specified point, given a */ /* discrete set of values of the function. */ /* Users of this routine must choose the number of points to use */ /* in their interpolation method. The authors of Reference [1] have */ /* this to say on the topic: */ /* Unless there is solid evidence that the interpolating function */ /* is close in form to the true function f, it is a good idea to */ /* be cautious about high-order interpolation. We */ /* enthusiastically endorse interpolations with 3 or 4 points, we */ /* are perhaps tolerant of 5 or 6; but we rarely go higher than */ /* that unless there is quite rigorous monitoring of estimated */ /* errors. */ /* The same authors offer this warning on the use of the */ /* interpolating function for extrapolation: */ /* ...the dangers of extrapolation cannot be overemphasized: */ /* An interpolating function, which is perforce an extrapolating */ /* function, will typically go berserk when the argument x is */ /* outside the range of tabulated values by more than the typical */ /* spacing of tabulated points. */ /* For Lagrange interpolation on unequally spaced abscissa values, */ /* see the SPICELIB routine LGRINT. */ /* $ Examples */ /* 1) Fit a cubic polynomial through the points */ /* ( -1, -2 ) */ /* ( 1, -8 ) */ /* ( 3, 26 ) */ /* ( 5, 148 ) */ /* and evaluate this polynomial at x = 2. */ /* PROGRAM TEST_LGRESP */ /* DOUBLE PRECISION LGRESP */ /* DOUBLE PRECISION ANSWER */ /* DOUBLE PRECISION FIRST */ /* DOUBLE PRECISION STEP */ /* DOUBLE PRECISION YVALS (4) */ /* DOUBLE PRECISION WORK (4) */ /* INTEGER N */ /* N = 4 */ /* FIRST = -1.D0 */ /* STEP = 2.D0 */ /* YVALS(1) = -2.D0 */ /* YVALS(2) = -8.D0 */ /* YVALS(3) = 26.D0 */ /* YVALS(4) = 148.D0 */ /* ANSWER = LGRESP ( N, FIRST, STEP, YVALS, WORK, 2.D0 ) */ /* WRITE (*,*) 'ANSWER = ', ANSWER */ /* END */ /* The returned value of ANSWER should be 1.D0, since the */ /* unique cubic polynomial that fits these points is */ /* 3 2 */ /* f(x) = x + 2x - 4x - 7 */ /* We also could have invoked LGRESP with the reference */ /* ANSWER = LGRESP ( N, FIRST, STEP, YVALS, YVALS, 2.D0 ) */ /* if we wished to; in this case YVALS would have been */ /* modified on output. */ /* If we had solved the same problem using a negative step, */ /* we would have set the elements of YVALS in reverse order, */ /* as shown below: */ /* FIRST = 5.D0 */ /* STEP = -2.D0 */ /* YVALS(1) = 148.D0 */ /* YVALS(2) = 26.D0 */ /* YVALS(3) = -8.D0 */ /* YVALS(4) = -2.D0 */ /* ANSWER = LGRESP ( N, FIRST, STEP, YVALS, WORK, 2.D0 ) */ /* The returned value of ANSWER would still be 1.D0. */ /* $ Restrictions */ /* None. */ /* $ Literature_References */ /* [1] "Numerical Recipes---The Art of Scientific Computing" by */ /* William H. Press, Brian P. Flannery, Saul A. Teukolsky, */ /* William T. Vetterling (see sections 3.0 and 3.1). */ /* $ Author_and_Institution */ /* N.J. Bachman (JPL) */ /* $ Version */ /* - SPICELIB Version 1.0.0, 14-AUG-1993 (NJB) */ /* -& */ /* $ Index_Entries */ /* interpolate function using Lagrange polynomial */ /* Lagrange interpolation */ /* -& */ /* SPICELIB functions */ /* Local variables */ /* Check in only if an error is detected. */ if (return_()) { ret_val = 0.; return ret_val; } /* No data, no interpolation. */ if (*n < 1) { ret_val = 0.; chkin_("LGRESP", (ftnlen)6); setmsg_("Array size must be positive; was #.", (ftnlen)35); errint_("#", n, (ftnlen)1); sigerr_("SPICE(INVALIDSIZE)", (ftnlen)18); chkout_("LGRESP", (ftnlen)6); return ret_val; } /* The step size must be non-zero. */ if (*step == 0.) { ret_val = 0.; chkin_("LGRESP", (ftnlen)6); setmsg_("Step size was zero.", (ftnlen)19); sigerr_("SPICE(INVALIDSTEPSIZE)", (ftnlen)22); chkout_("LGRESP", (ftnlen)6); return ret_val; } /* We can simplify the interpolation problem by shifting */ /* and scaling the abscissa values so that they start at 1 */ /* and are separated by a unit step. All we need to do is */ /* shift and scale X. */ newx = (*x - *first) / *step + 1.; /* We're going to compute the value of our interpolating polynomial */ /* at X by taking advantage of a recursion relation between */ /* Lagrange polynomials of order n+1 and order n. The method works */ /* as follows for the case of abscissa values that are not */ /* necessarily uniformly spaced: */ /* Define */ /* P (x) */ /* i(i+1)...(i+j) */ /* to be the unique Lagrange polynomial that interpolates our */ /* input function at the abscissa values */ /* x , x , ... x . */ /* i i+1 i+j */ /* Then we have the recursion relation */ /* P (x) = */ /* i(i+1)...(i+j) */ /* x - x */ /* i */ /* ----------- * P (x) */ /* x - x (i+1)...(i+j) */ /* i i+j */ /* x - x */ /* i+j */ /* + ----------- * P (x) */ /* x - x i(i+1)...(i+j-1) */ /* i i+j */ /* Repeated application of this relation allows us to build */ /* successive columns, in left-to-right order, of the */ /* triangular table */ /* P (x) */ /* 1 */ /* P (x) */ /* 12 */ /* P (x) P (x) */ /* 2 123 */ /* P (x) */ /* 23 . */ /* P (x) */ /* . 234 . */ /* . */ /* . . . */ /* . */ /* . . P (x) */ /* . . 12...N */ /* . */ /* . */ /* . */ /* P (x) */ /* (N-2)(N-1)N */ /* P (x) */ /* (N-1)N */ /* P (x) */ /* N */ /* and after N-1 steps arrive at our desired result, */ /* P (x). */ /* 12...N */ /* In the current case, we've arranged the problem so that */ /* x = i. */ /* i */ /* We'll use the scratch array WORK to contain the current column of */ /* our interpolation table. To start out with, WORK(I) will contain */ /* P (x). */ /* I */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { work[i__ - 1] = yvals[i__ - 1]; } /* Compute columns 2 through N of the table. */ i__1 = *n - 1; for (j = 1; j <= i__1; ++j) { i__2 = *n - j; for (i__ = 1; i__ <= i__2; ++i__) { c1 = (doublereal) (i__ + j) - newx; c2 = newx - (doublereal) i__; work[i__ - 1] = (c1 * work[i__ - 1] + c2 * work[i__]) / ( doublereal) j; } } /* Our result is sitting in WORK(1) at this point. */ ret_val = work[0]; return ret_val; } /* lgresp_ */