SUBROUTINE CSPLN(X,Y,C,N)
*     Calculation of spline coefficients from  N pairs x,y in 
*     arrays X(N) and Y(N) (array X in increasing order), using the
*     not-a-knot end conditions.
*     The array C, dimensioned C(5,N), contains after return
*        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)
*     and can be used for interpolation.
*     Calculation of interpolated value for given X 
*        with X(I) =< X =< X(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,*),X(*),Y(*)
      IF(N.LT.3) RETURN
      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
*     not-a-knot at left side
      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 
      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
*     not-a-knot at right side
      IF(N.EQ.3) THEN
         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 only)
      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)
      END