Cubic Spline Curve Fitting

Oforth examples. Feel free to post your own code.

Cubic Spline Curve Fitting

Postby Doug » 20 Nov 2018 16:16

I have found the following code to be very useful for highly accurate curve fitting. In particular I have used the algorithm to create airfoil templates for full size aircraft and for the keel of a 10 meter sailboat.

The algorithm is based on the following old unstructured Fortran code. After that I have included the routine converted to Oforth.

SUBROUTINE AKIMA ( X, Y, N, XX, YY, SLP )
C
C
C...PURPOSE General-purpose monovariate interpolation routine
C using a locally-fitted cubic spline. One point is
C interpolated from the input data arrays.
C
C...INPUT: X An input array of abscissas in ascending or
C descending value order.
C Y An input array of corresponding ordinates. The
C function Y(X) must be single-valued.
C N Size for the 'X' and 'Y' arrays. 'N' is an integer.
C (the input arrays may be any size)
C
C XX Input 'X' point at which an interpolated 'Y' value is
C desired.
C
C...OUTPUT: YY Interpolated 'Y' ordinate.
C SLP Interpolated slope (DY/DX).
C
C...DISCUSSION This spline method produces an interpolated curve that
C is relatively free of the oscillatory behavior normally
C associated with cubic splines. The curve produced will
C be continuous in 'Y' and 'DY/DX' but further derivatives
C may be discontinuous. The interpolated slope should be
C treated with some caution, however, due to the tendency
C of this spline curve to concentrate changes of curvature
C at the input points (interval ends). If only 2 points are
C input the curve will be linear, if 3 points are given
C the curve will be parabolic, more than 3 points will
C produce a cubic. Extrapolation beyond the bounds of the
C input data is done with a quadratic through the last 3
C points at that data boundary.
C
C This routine is intended as a replacement for B.Wainfan's
C AKIMAD. It is shorter, twice as fast and it works on
C input data in ascending or descending order. The calling
C sequence is more natural and is easier to use for the
C bulk of applications.
C
C The coding is compatible with IBM or VAX Fortran 77
C and IBM Fortran IV G or H with optimization.
C
C...ORIGIN Harold Youngren CALAC Dept. 72-71 3/81
C
C...REFERENCE Akima, Hiroshi, "A New Method of Interpolation and
C Smooth Curve Fitting Based on Local Procedures",
C Journal of the Association for Computing Machines,
C Vol. 17, No. 4, Oct 1970, pages 589-602
C
C Wainfan, B. S., "The Akima Subroutines: Nonlinear
C Interpolation by Local Polynomial Fit", LR 29244,
C Oct 30, 1979.
C
C
DIMENSION X(N), Y(N), D(5), T(2)
C
C
C...Check for a degenerate case ( X(1)=X(N) ).
C
IF (X(1) .NE. X(N)) GO TO 10
YY = Y(1)
SLP = 0.
GO TO 70
C
C
C...Find which interval contains the point by binary search.
C...The binary search loop is terminated when the search residual
C...(NSTEP) is zero. The index 'I' will point to the input point
C...lower than or equal to the desired point for the 'X' values
C...in ascending order or to the input point greater than or equal
C...to the desired point for descending order 'X' values.
C
10 XORDR = 1.0
IF (X(1) .GT. X(N)) XORDR = -1.0
C
IBOT = 1
ITOP = N
XXO = XX * XORDR
C
20 NSTEP = ( ITOP - IBOT ) / 2
I = IBOT + NSTEP
XO = X(I) * XORDR
IF ( XXO .GE. XO ) IBOT = I
IF ( XXO .LT. XO ) ITOP = I
IF ( NSTEP .NE. 0 ) GO TO 20
C
C
C...Calculate the straight line slopes between adjacent input points.
C...D(3) is the slope on the interval of interpolation. If the other
C...slopes D(1), D(2), D(4) OR D(5) are not defined they will be
C...created by quadratic extrapolation (only at start and end of data).
C
DO 30 J = 1, 5
K = I + (J-2)
IF ( ((K-1) .GE. 1) .AND. (K .LE. N) )
* D(J) = ( Y(K) - Y(K-1) ) / ( X(K) - X(K-1) )
30 CONTINUE
C
C...Synthesize upper and lower slopes if required. Check for
C...single line segment input (N=2).
C
IF (N .EQ. 2) D(2) = D(3)
C
IF ((I+2) .GT. N) D(4) = 2. * D(3) - D(2)
IF ((I+3) .GT. N) D(5) = 2. * D(4) - D(3)
IF ((I-1) .LT. 1) D(2) = 2. * D(3) - D(4)
IF ((I-2) .LT. 1) D(1) = 2. * D(2) - D(3)
C
C
C...Calculate the slopes (T(1),T(2)) at the lower and upper
C...points bounding the interval of interpolation. If the point is
C...at an intersection of straight line segments the slope is
C...defined by the average of the adjacent segment slopes.
C
DO 50 J = 1, 2
A = ABS( D(J+3) - D(J+2) )
B = ABS( D(J+1) - D(J) )
IF ((A + B) .NE. 0.) GO TO 40
A = 1.
B = 1.
40 T(J) = ( A*D(J+1) + B*D(J+2) ) / ( A + B )
50 CONTINUE
C
C
C...Check if desired point is on upper point of interval. This
C...reduces error at the transition points between intervals.
C
IF (XX .NE. X(I+1)) GO TO 60
YY = Y(I+1)
SLP = T(2)
GO TO 70
C
C...Calculate the cubic coefficients.
C
60 XINT = X(I+1) - X(I)
XDIF = XX - X(I)
P0 = Y(I)
P1 = T(1)
P2 = ( 3.*D(3) - 2.*T(1) - T(2) ) / XINT
P3 = ( T(1) + T(2) - 2.*D(3) ) / (XINT*XINT)
C
C...Calculate the Y-value and the slope.
C
YY = P0 + XDIF*( P1 + XDIF*( P2 + XDIF*P3 ) )
SLP = P1 + XDIF*( 2.*P2 + XDIF*( 3.*P3 ) )
C
70 RETURN
END

Code: Select all


\ Airfoil plotting is one use of the Akima curve fitting program.
\ A very complete database of airfoil data for downloading can be found
\ at the following site:
\ https://m-selig.ae.illinois.edu/ads/coord_database.html



tvar: x \ container for an array 
tvar: y \ container for an array

 
\ Example data, NACA 0010-64 top surface only (it is a symmetric airfoil)
\ It would not be difficult to add Oforth code to read the x,y points directly
\ from the downloaded data files.

[ 1.0000, 0.9500, 0.9000, 0.8000, 0.7000, 0.6000, 0.5000, 0.4000, 0.3000, 0.2000, 0.1500, 0.1000, 0.0750, 0.0500, 0.0250, 0.0125, 0.0000 ]
to x
[ 0.00100, 0.00856, 0.01556, 0.02767, 0.03733, 0.04433, 0.04856, 0.05000, 0.04856, 0.04411, 0.04056, 0.03533, 0.03178, 0.02722, 0.02044, 0.01511, 0.00000 ]
to y
 
 
\ ...Check for a degenerate case ( X(1)=X(N) ).
\ ...and for equal size X, Y arrays
\ ...and if less than 5 points

: check ( n -- b )
  n 5 < if "Less than 5 points" abort then
  y size n <> if "arrays x and y not same size" abort then
  1 x at  n x at = ;

\ ...Find which interval contains the point by binary search.
\ ...The binary search loop is terminated when the search residual
\ ...(NSTEP) is zero. The index 'I' will point to the input point
\ ...lower than or equal to the desired point for the 'X' values
\ ...in ascending order or to the input point greater than or equal
\ ...to the desired point for descending order 'X' values.

: interval ( xx i n --  i )
  | xord xxo xo itop ibot nstep |
  1 -> xord
  1 x at  n x at > if -1 -> xord then
 
  1 -> ibot
  n -> itop
  xx xord * -> xxo

      begin
       itop ibot - 2 / -> nstep
       ibot nstep + -> i
       xord i x at * -> xo
       xxo xo >= if i -> ibot else i -> itop then
       nstep 0 = 
     until i
;

\ ...Calculate the straight line slopes between adjacent input points.
\ ...D(3) is the slope on the interval of interpolation. If the other
\ ...slopes D(1), D(2), D(4) OR D(5) are not defined they will be
\ ...created by quadratic extrapolation (only at start and end of data).

: slopes ( i n d -- d ) \ it may be possible to skip the step of returning the d object to the calling function
  | j k |
  1 5 for: j [ j 2 - i + -> k
               k 1 - 1 >=  k n <=  and
               if k y at  k 1 - y at -  k x at  k 1 - x at  - / j swap d put then ]

\ ...Synthesize upper and lower slopes if required.
       
        i 2 + n > if 2 3 d at *  2 d at -  4 swap d put then
        i 3 + n > if 2 4 d at *  3 d at -  5 swap d put then
        i 1 - 1 < if 2 3 d at *  4 d at -  2 swap d put then
        i 2 - 1 < if 2 2 d at *  3 d at -  1 swap d put then
        d
;


\ INPUT ARRAYS     X  An input array of at least 5 abscissas in ascending or
\                     descending value order.
\                  Y  An input array of corresponding ordinates. The
\                     function Y(X) must be single-valued.
 
\ ...INPUT:    XX  Input 'X' point at which an interpolated 'Y' value is
\                  desired.
\
\ ...OUTPUT:   YY  Interpolated 'Y' ordinate.
\             SLP  Interpolated slope at point 'XX,YY'(DY/DX).

: akima ( xx -- yy slp )
  | i j n a b d t xint xdif p0 p1 p2 p3 |
  2 Array newSize -> t
  5 Array newSize -> d
  x size -> n
  n check if 1 t at 0.0 return then
  xx i n interval -> i
  i n d slopes -> d 

\ ...Calculate the slopes (T(1),T(2)) at the lower and upper
\ ...points bounding the interval of interpolation. If the point is
\ ...at an intersection of straight line segments the slope is
\ ...defined by the average of the adjacent segment slopes.

  1 2 for: j [ j 3 + d at  j 2 + d at - abs -> a
               j 1 + d at  j d at - abs -> b
               a b + 0 =  if
                            1 -> a
                            1 -> b
                          then
               a j 1 + d at  *  b j 2 + d at * + a b + /  j swap t put ]


\ ...Check if desired point is on upper point of interval. This
\ ...reduces error at the transition points between intervals.

        xx i 1 + x at <>
        if
\ ...Calculate the cubic coefficients.
           i 1 + x at  i x at - -> xint
           xx i x at - -> xdif
           i y at -> p0
           1 t at -> p1
           3 3 d at *  2 1 t at * -  2 t at - xint / -> p2
           1 t at 2 t at +  2 3 d at * -  xint xint * / -> p3
                     
\ ...Calculate the Y-value and the slope.
           xdif p3 * p2 +  xdif *  p1 + xdif * p0 + \ yy
           3 p3 * xdif * 2 p2 * + xdif * p1 +  \ slp
           
        else
           i 1 + y at \ yy
           2 t at \ slp
        then
 ;

Doug
 
Posts: 27
Joined: 20 Jul 2018 14:26

Return to Oforth examples

Who is online

Users browsing this forum: No registered users and 0 guests

cron