SUBROUTINE BSKFIT(X,Y,W,N,K,T,NKN,C)
*     Least squares spline fit of B-spline of order K to 
*     data Y in array Y(N), with abscissa values X in arrays X(N)
*     and weights W in array W(N). 
*     The array T(NKN) contains NKN (ordered) knot positions.
*     The resulting B-spline coefficients are in an array of dimension 
*     NDIM = NKN-K, i.e. C(NKN-K).
*     Limitation: NDIM <= ND=100 ; K <= 10 
      REAL X(N),Y(N),W(N), T(NKN),C(*)   ! C(NKN-K)
      PARAMETER (ND=100)
      REAL Q(4*ND),B(10)
*     ...
      DO I=1,NKN-K                      ! reset normal equation
       C(I)=0.0 
      END DO 
      DO I=1,K*(NKN-K)
       Q(I)=0.0
      END DO 
      DO M=1,N
       L=MAX(K,MIN(LEFT(X(M),T,NKN),NKN-K)) !index within limits
       CALL BSPLNK(X(M),K,B,L,T)       ! calculate B-splines
       DO I=1,K                        ! add to normal equation
        C(I+L-K)=C(I+L-K)+B(I)*Y(M)*W(M)
        IQ=1+K*(I+L-K)-I-K
        DO J=I,K
         Q(J+IQ)=Q(J+IQ)+B(I)*B(J)*W(M)    
        END DO
       END DO 
      END DO 
      CALL BMSOL(Q,C,K,NKN-K)       ! solve band matrix equation
      END