
      SUBROUTINE CSPNK(X,Y,C,N)
*
*     Calculation of spline coefficients from pairs x,y
*        with not-a-knot end-conditions
*
*     Spline coefficients C(1...5,N) of cubic spline F(X) are:
*
*          C(1,I) =  F(X(I)) = Y(I)
*          C(2,I) =  first derivative  F'(X(I))
*          C(3,I) =  second derivative F''(X(I))/2.0
*          C(4,I) =  third derivative  F'''(X(I))/6.0
*          C(5,I) =  X(I)                          for I=1...N
*
*     Calculation of interpolated value for given X
*        with C(5,I)=X(I) =< X =< X(I+1)=C(5,I+1):
*          H=X-C(5,I)
*          F=C(1,I)+H*(C(2,I)+H*(C(3,I)+H*C(4,I)))
*     ------------------------------------------------------------------
*
      REAL C(5,N),X(*),Y(*)
*     ...
      IF(N.LT.3) GOTO 100
      DO I=1,N
       C(1,I)=Y(I)
       C(5,I)=X(I)
      END DO
*     first differences of x and first divided differences of data
      DO I=2,N
       C(3,I)= C(5,I)-C(5,I-1)
       C(4,I)=(C(1,I)-C(1,I-1))/C(3,I)
      END DO
      C(4,1)=C(3,3)
      C(3,1)=C(3,2)+C(3,3)
      C(2,1)=((C(3,2)+2.0*C(3,1))*C(4,2)*C(3,3)+C(3,2)**2*C(4,3))/C(3,1)
*     generate equations and carry out forward pass of Gauss elimination
      DO I=2,N-1
       G=-C(3,I+1)/C(4,I-1)
       C(2,I)=G*C(2,I-1)+3.0*(C(3,I)*C(4,I+1)+C(3,I+1)*C(4,I))
       C(4,I)=G*C(3,I-1)+2.0*(C(3,I)+C(3,I+1))
      END DO
      IF(N.EQ.3) THEN
*        special case of three points
         C(2,N)=2.0*C(4,N)
         C(4,N)=1.0
         G=-1.0/C(4,N-1)
      ELSE
         G=C(3,N-1)+C(3,N)
         C(2,N)=((C(3,N)+2.0*G)*C(4,N)*C(3,N-1)+C(3,N)**2*(C(1,N-1)
     +          -C(1,N-2))/C(3,N-1))/G
         G=-G/C(4,N-1)
         C(4,N)=C(3,N-1)
      END IF
*     complete forward pass
      C(4,N)=C(4,N)+G*C(3,N-1)
      C(2,N)=(G*C(2,N-1)+C(2,N))/C(4,N)
*     backward substitution
      DO I=N-1,1,-1
       C(2,I)=(C(2,I)-C(3,I)*C(2,I+1))/C(4,I)
      END DO
*     generate coefficients
      DO I=2,N
       DX=C(3,I)
       DVD1=(C(1,I)-C(1,I-1))/DX
       DVD3=C(2,I-1)+C(2,I)-2.0*DVD1
       C(3,I-1)=(DVD1-C(2,I-1)-DVD3)/DX
       C(4,I-1)=(DVD3/DX)/DX
      END DO
*     calculation of coefficients at X(N) (for completeness)
      DX=C(5,N)-C(5,N-1)
      C(2,N)=C(2,N-1)+DX*(2.0*C(3,N-1)+DX*3.0*C(4,N-1))
      C(3,N)=C(3,N-1)+DX*3.0*C(4,N-1)
      C(4,N)=C(4,N-1)
  100 RETURN
      END

      FUNCTION FCSPNK(X,C,N)
*     value of spline function C(5,N) for given value of X
      REAL C(5,N)
*     ...
      DO I=1,N-1
       IF(X.LT.C(5,I)) GOTO 10
      END DO
      I=N+1
   10 I=MAX(1,I-1)
      H=X-C(5,I)
      FCSPNK=C(1,I)+H*(C(2,I)+H*(C(3,I)+H*C(4,I)))
      RETURN
      END
