stellinstall/nspline_8f_source.html

C******************** START FILE SPLINE.FOR ; GROUP TRKRLIB ******************

CC

CCCCCCCCCCCCCCCC

C   THE CODES (SPLINE & SEVAL) ARE TAKEN FROM:

C   FORSYTHE,MALCOLM AND MOLER, "COMPUTER METHODS FOR

C   MATHEMATICAL COMPUTATIONS",PRENTICE-HALL, 1977.

C

C   THE CODES (SPLEEN,SPLAAN & SPEVAL) ARE ADAPTATIONS

C   BY R.M. WIELAND FOR SPECIAL CASES ... SEE COMMENTS

c

c  this is nspline -- natural cubic spline

c    second derivative of interpolant set to zero at endpts

c


      SUBROUTINE nspline (N, X, Y, B, C, D)

      INTEGER N

      REAL X(N), Y(N), B(N), C(N), D(N)

C

C  THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED

C  FOR A CUBIC INTERPOLATING SPLINE

C

C    S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3

C

C    FOR  X(I) .LE. X .LE. X(I+1)

C

C  INPUT..

C

C    N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)

C    X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER

C    Y = THE ORDINATES OF THE KNOTS

C

C  OUTPUT..

C

C    B, C, D  = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.

C

C  USING  P  TO DENOTE DIFFERENTIATION,

C

C    Y(I) = S(X(I))

C    B(I) = SP(X(I))

C    C(I) = SPP(X(I))/2

C    D(I) = SPPP(X(I))/6  (DERIVATIVE FROM THE RIGHT)

C

CCCCCCCCCCCCCCC

C  THE ACCOMPANYING FUNCTION SUBPROGRAM  SEVAL  CAN BE USED

C  TO EVALUATE THE SPLINE.

C

C

      INTEGER NM1, NM2, IB, I

      REAL T

C

      nm1 = n-1

      IF ( n .LT. 2 ) RETURN

C

C

C DMC - CHECK ORDINATES

      DO 5 i=1,nm1

        ip1=i+1

        IF(x(ip1).LE.x(i)) then

             write(6,9000)

             write(6,9001) (x(ii),ii=1,n)

             zbomb=sqrt(-1.0/(x(1)*x(1)))

             write(6,*) zbomb

             write(6,*) ' nspline -- code stop.'

             call exit(1)

 9000        FORMAT(/

     >' ? UNORDERED ORDINATE ARRAY PASSED TO SPLINE SUBROUTINE:')

 9001        FORMAT(1x,5(1x,1pe12.5))

        endif

 5    CONTINUE

C

      IF ( n .LT. 3 ) GO TO 50

C

C  SET UP TRIDIAGONAL SYSTEM

C

C  B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE.

C

      d(1) = x(2) - x(1)

      c(2) = (y(2) - y(1))/d(1)

      DO 10 i = 2, nm1

         d(i) = x(i+1) - x(i)

         b(i) = 2.*(d(i-1) + d(i))

         c(i+1) = (y(i+1) - y(i))/d(i)

         c(i) = c(i+1) - c(i)

   10 CONTINUE

C

C  END CONDITIONS:  2nd deriv = 0 at endpts; solve reduced n-2 x n-2

C  system.

      c(1) = 0.

      c(n) = 0.

C

C  FORWARD ELIMINATION

C

   15 DO 20 i = 3, nm1

         t = d(i-1)/b(i-1)

         b(i) = b(i) - t*d(i-1)

         c(i) = c(i) - t*c(i-1)

   20 CONTINUE

C

C  BACK SUBSTITUTION

C

      nm2=nm1-1

      c(nm1) = c(nm1)/b(nm1)

      DO 30 ib = 2, nm2

         i = n-ib

         c(i) = (c(i) - d(i)*c(i+1))/b(i)

   30 CONTINUE

C

C  C(I) IS NOW THE SIGMA(I) OF THE TEXT

C

C  COMPUTE POLYNOMIAL COEFFICIENTS

C

      b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2.*c(n))

      DO 40 i = 1, nm1

         b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.*c(i))

         d(i) = (c(i+1) - c(i))/d(i)

         c(i) = 3.*c(i)

   40 CONTINUE

      c(n) = 3.*c(n)

      d(n) = d(n-1)

      RETURN

C

   50 b(1) = (y(2)-y(1))/(x(2)-x(1))

      c(1) = 0.

      d(1) = 0.

      b(2) = b(1)

      c(2) = 0.

      d(2) = 0.

      RETURN

      END

C******************** END FILE SPLINE.FOR ; GROUP TRKRLIB ******************