stellinstall/pspline_8f_source.html

C  PSPLINE -- modified from SPLINE.FOR -- dmc 3 Oct 1995

C

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

C   FORSYTHE,MALCOLM AND MOLER, "COMPUTER METHODS FOR

C   MATHEMATICAL COMPUTATIONS",PRENTICE-HALL, 1977.

c

c  PSPLINE -- periodic spline -- adaptation of SPLINE.FOR by D. McCune

c  October 1995:  The function to be splined is taken as having periodic

C  derivatives (it need not be periodic itself), so that

c  the interpolating function satisfies

c            y'(x(1))=y'(x(N))

c           y''(x(1))=y''(x(N))

c  this results in a fully specified system (no addl BC assumptions

c  required); it is not a tridiagonal system but a modified NM1xNM1

c  tridiagonal system with non-zero offdiagonal corner (1,NM1), (NM1,1)

c  elements.  It remains symmetric & diagonally dominant, and is

c  solved by Gaussian elimination w/o pivoting.  NM1=N-1.

c

      SUBROUTINE pspline (N, X, Y, B, C, D, WK)

      INTEGER N

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

C

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

C  FOR A periodic 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    WK(...) is used as a workspace during the Guassian elimination.

C    it represents the rightmost column of the array

C

C    note B(N),C(N),D(N) give coeffs for a periodic continuation into

C    the next period, up to X(N)+(X(2)-X(1)) only.

C

C    SEVAL can be used to evaluate the spline, but, it is up to the

C    SEVAL caller to bring the argument into the normal range X(1) to X(N).

C

C  INPUT..

C

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

C        includes a complete period of the data, the last point being

C        a repeat of the first point.

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

C        X(N)-X(1) is the period of the periodic function represented.

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

C-------------------------------------

C

      nm1 = n-1

      nm2 = nm1-1

C

      IF ( n .LT. 4 ) THEN

         write(6,9901)

 9901    format(/

     >   ' ?PSPLINE -- at least 4 pts required for periodic spline.')

         return

      ENDIF

C

C  SET UP MODIFIED NM1 x NM1 TRIDIAGONAL SYSTEM:

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

C  WK(1:NM2) = rightmost column above diagonal

C  WK(NM1)   = lower left corner element

C

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

      d(nm1) = x(n) - x(nm1)

      b(1) = 2.0*(d(1)+d(nm1))

      c(1) = (y(n) - y(nm1))/d(nm1)

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

      c(1) = c(2) - c(1)

      wk(1) = d(nm1)

      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)

         wk(i) = 0.0

   10 CONTINUE

      wk(nm2) = d(nm2)

      wk(nm1) = d(nm1)

C

C  END CONDITIONS -- implied by periodicity

C     C(1)=C(N)  B(1)=B(N)  D(1)=D(N) -- no additional assumption needed.

C

C  FORWARD ELIMINATION

C   WK(1)--WK(NM2) represent the rightmost column above the

C   diagonal; WK(NM1) represents the non-zero lower left corner element

C   which in each step is moved one column to the right.

C

   15 DO 20 i = 2, nm2

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

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

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

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

         q = wk(nm1)/b(i-1)

         wk(nm1) = -q*d(i-1)

         b(nm1) = b(nm1) - q*wk(i-1)

         c(nm1) = c(nm1) - q*c(i-1)

   20 CONTINUE

C

C  correct the (NM1,NM2) element

C

      wk(nm1) = wk(nm1) + d(nm2)

C

C  complete the forward elimination:  now WK(NM1) and WK(NM2) are

C  the off diagonal elements of the 2x2 at the lower right hand corner

C

      t = wk(nm1)/b(nm2)

      b(nm1) = b(nm1) - t*wk(nm2)

      c(nm1) = c(nm1) - t*c(nm2)

C

C  BACK SUBSTITUTION

C

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

      c(nm2) = (c(nm2) - wk(nm2)*c(nm1))/b(nm2)

C

      DO 30 ib = 3, nm1

         i = n-ib

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

   30 CONTINUE

C

C  PERIODICITY:

C

      c(n)=c(1)

C

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

C

C  COMPUTE POLYNOMIAL COEFFICIENTS

C

      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 periodic continuation...

      b(n) = b(1)

      c(n) = c(1)

      d(n) = d(1)

C

      RETURN

      END