REAL FUNCTION POLY( X, XPTS, YPTS, NDEG )
INTEGER, INTENT(IN ) :: NDEG ! degree to use
REAL , INTENT(IN ) :: X ! X-value to compute Y for
REAL , INTENT(IN ) :: XPTS ( NDEG + 1 ) ! table of X-values to be used
REAL , INTENT(IN ) :: YPTS ( NDEG + 1 ) ! table of Y-values to interpolate from
float POLY( const float *x,
const float *xpts,
const float *ypts,
const int *n ) ;
X on the curve determined by
XPTS and YPTS using Newton
divided-differences. Although strictly speaking no ordering
assumptions are required, numerical stability is greatest
(round-off error is minimized) if the XPTS are in
sorted order, and X lies near the middle of the list.
When interpolating from a large table, it is the
responsibility of the caller to pass the portion
of the table appropriate for the X at which the
interpolation is being evaluated. Note that high-order polynomial
interpolation is subject to Gibbs-Phenomenon numerical errors, and
is to be avoided.
For Fortran-90 declarations and interface checking:
USE M3UTILIO
NOTE: high-order purely polynomial interpolations have stability problems. NDEG <= 5 is recommended. -- CJC
See also PCOEF().
#include "iodecl3.h" if called from C.
X should lie in the domain specified by the array
XPTS of X-values so that
POLY() is doing interpolation instead of extrapolation.
X is between XPT(17) and
XPT(18), and you want to
do a degree-3 interpolation (requiring XPT(16) ... XPT(19) ):
...
USE M3UTILIO
...
REAL X, Y
REAL XPT( 100 )
DATA XPT / ... /
REAL YPT( 100 )
DATA YPT / ... /
...
Y = POLY( X, XPT( 16 ), YPT( 16 ), 3 ) ! interpolate Y-values to X
...
...
#include "iodecl3.h"
...
float x, y, xpt[], ypt[] ;
int deg ;
...
y = POLY( &x, &xpt, &ypt, ° ) ;
/* Now y is result of degree-deg polynomial interpolation of curve
determined by xpt[] and ypt[] to x */
...
To: Models-3/EDSS I/O API: The Help Pages