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

      FUNCTION CSPF(X,C,N)
*     Calculate function value for given argument X, for spline
*     function in array C, dimensioned C(5,N), obtained by previous
*     call of CSPLN.
      REAL C(5,N)
      DATA  I/1/
      IF(I.GE.N) I=(N+1)/2
      IF(X.GE.C(5,I)) THEN   ! make use of previous index 
         IH=MIN(I+2,N-1)
         IF(X.GE.C(5,IH)) IH=N-1         
         IL=I
      ELSE 
         IL=MAX(I-2,1)
         IF(X.LT.C(5,IL)) IL=1
         IH=I
      END IF 
   10 I=(IL+IH)/2     ! binary search 
      IF(I.NE.IL) THEN
         IF(X.GE.C(5,I)) THEN
            IL=I
         ELSE
            IH=I
         END IF
         GOTO 10
      END IF
      H=X-C(5,I)
      CSPF=C(1,I)+H*(C(2,I)+H*(C(3,I)+H*C(4,I)))
      END

      SUBROUTINE USPCN(Y,A,N,G,H)
*     Calculation of B-spline coefficients A(N) from equidistant function
*     values in array Y(N), using the not-a-knot condition.
*     Arrays Y(*) and A(*) may be identical. 
*     G(N) and H(N) are auxiliary arrays used during computation. 
      REAL Y(N),A(N),G(N),H(N)
      DO I=2,N
       G(I)=Y(I)-Y(I-1)
      END DO
      H(1)=(5.0*G(2)+G(3))/2.0
      H(2)=-H(1)+3.0*(G(3)+G(2))
      G(2)=2.0
      DO I=3,N-1
       H(I)=-H(I-1)/G(I-1)+3.0*(G(I+1)+G(I))
       G(I)=4.0-1.0/G(I-1)
      END DO
      H(N)=(5.0*G(N)+Y(N-1)-Y(N-2))/2.0
      G(N)=1.0-2.0/G(N-1)
*     backward ...
      H(N)=(-2.0*H(N-1)/G(N-1)+H(N))/G(N)
      DO I=N-1,2,-1
       H(I)=(H(I)-H(I+1))/G(I)
      END DO
      H(1)=(H(1)-2.0*H(2))
*     ... and forward solution
      ALN=Y(N)+Y(N)-Y(N-1)-(2.0*H(N)+H(N-1))/3.0
      DO I=1,N-1
       A(I)=Y(I)+Y(I)-Y(I+1)+(2.0*H(I)+H(I+1))/3.0
      END DO
      A(N)=ALN
      END 

      FUNCTION USPF(X,A,N,XA,XB)
*     Returned is the function value at X of B-spline function.
*     A(N) is array of coefficients, obtained from USPCN.
*     XA and  XB are X-values of first and last function value 
*     Extrapolation is done for X outside the range XA ... XB.
      REAL A(N)
      R=FLOAT(N-1)*(X-XA)/(XB-XA)
      I=MAX(1,MIN(INT(R),N-3))
      Z=R-FLOAT(I)
      USPF
     +    =(A(I)+4.0*A(I+1)+A(I+2)+Z*(3.0*(-A(I)+A(I+2))
     +     +Z*(3.0*(A(I)-2.0*A(I+1)+A(I+2))
     +     +Z*(-A(I)+3.0*A(I+1)-3.0*A(I+2)+A(I+3)))))/6.0
      END

      SUBROUTINE BSPLNK(X,K,B,L,T)
*     Returned are B(1)...B(K) of B-spline values for given X and order K 
*     K = order of spline (1=const, 2=linear...) up to limit K=KMAX 
*     The array T(*) contains the knot positions.
*     In T(K-1) =<  x  < T(M+2-K) the K B-spline values add up to 1. 
*     Multiple knots allowed: (K - mult)-th derivative is discontinuous.
      PARAMETER (KMAX=10)
      REAL B(K),DL(KMAX-1),DR(KMAX-1),SV,TR,T(*)
      IF(K.GT.KMAX) STOP ' BSPLNK: K too large '
      B(1)=1.0
      DO J=1,K-1
       DR(J)=T(L+J)-X
       DL(J)=X-T(L-J+1)
       SV=0.0
       DO I=1,J
        TR  =B(I)/(DR(I)+DL(J-I+1))
        B(I)=SV+TR*DR(I)
        SV  =   TR*DL(J-I+1)
       END DO
       B(J+1)=SV
       IF(B(  1)+1.0.EQ.1.0) B(  1)=0.0 ! avoid underflows
       IF(B(J+1)+1.0.EQ.1.0) B(J+1)=0.0
      END DO
      END

      SUBROUTINE BSKFIT(X,Y,W,N,KS,T,NKN,C)
*     Least squares spline fit of B-spline of order KS 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 array of dimension 
*     NDIM = NKN-KS+1, i.e. C(NKN-KS+1).
*     Limitation: NDIM <= ND; KS <= 10 
      REAL X(N),Y(N),W(N), T(NKN),C(*)   ! C(NKN-KS+1)
      PARAMETER (ND=100)
      REAL Q(4*ND),B(10)
      DO I=1,NKN-KS+1                    ! reset normal equation
       C(I)=0.0 
      END DO 
      DO I=1,KS*(NKN-KS+1)
       Q(I)=0.0
      END DO 
      DO K=1,N
       L=MAX(KS-1,MIN(LEFT(X(K),T,NKN),NKN-KS+1)) !index within limits
       CALL BSPLNK(X(K),KS,B,L,T)       ! calculate B-splines
       DO I=1,KS                        ! add to normal equation
        C(I+L-KS+1)=C(I+L-KS+1)+B(I)*Y(K)*W(K)
        IQ=1+KS*(I+L-KS)-I
        DO J=I,KS
         Q(J+IQ)=Q(J+IQ)+B(I)*B(J)*W(K)    
        END DO
       END DO 
      END DO 
      CALL BMSOL(Q,C,KS,NKN-KS+1)       ! solve band matrix equation
      END

      FUNCTION BSKFUN(X,K,T,NKN,C)
*     Returned is the function value at X of B-spline function of 
*     order K within the knot range.
*     C  is the array of coefficients, and T(NKN) is the knot array.
      PARAMETER (KMAX=10) 
      REAL T(NKN),B(KMAX),C(*)
      DOUBLE PRECISION SUM
      L=MAX(K-1,MIN(LEFT(X,T,NKN),NKN-K+1)) ! index within limits
      CALL BSPLNK(X,K,B,L,T)                ! calculate four B-splines 
      SUM=B(1)*C(L-K+2)
      DO I=2,K
       SUM=SUM+B(I)*C(I+L-K+1) 
      END DO 
      BSKFUN=SUM
      END  

      SUBROUTINE ORTHFT(X,Y,W,N,NPAR,Q,P)
*     Construction of orthogonal polynomials from   
*     data Y in array Y(N), with abscissa values X in arrays X(N)
*     and weights W in array W(N) with NPAR parameters (up to order
*     NPAR-1 of the polynomial.
*     Q(2,N) is an auxiliary array, used during computation.
*     The array P, dimensioned P(5,NPAR), contains after return
*     the result of  the fit, to be used in subroutine ORTHFE later.
*        P(1,i) = alpha               P(1,1) = 100*N + NPAR
*        P(2,i) = beta                P(2,1) = variance estimate
*        P(3,i) = gamma               P(3,1) = p0(x) (constant)
*        P(4,i) = coefficient
*        P(5,i) = remaining chi**2
      REAL X(N),Y(N),W(N),Q(2,N),P(5,NPAR)
*     Double precision for the accumulation of sums
      DOUBLE PRECISION ABGDS(5),AL,BE,GA,DE,SQSUM
      EQUIVALENCE (ABGDS(1),AL),(ABGDS(2),BE),(ABGDS(3),GA),
     +            (ABGDS(4),DE),(ABGDS(5),SQSUM)
      NP=NPAR
      DO J=1,5
       ABGDS(J)=0.0D0
      END DO 
*     determination of constant term
      DO I=1,N
       GA=GA+W(I)       ! sum of weights
       DE=DE+W(I)*Y(I)  ! weighted y sum 
      END DO 
      IF(GA.NE.0.0) GA=1.0D0/SQRT(GA)
      DE  =GA*DE
      DO I=1,N
       SQSUM =SQSUM+W(I)*(Y(I)-GA*DE)**2 ! sum of squared residuals
       Q(1,I)=GA        ! initialize Q array
       Q(2,I)=0.0
      END DO
      DO J=1,5
       P(J,1)=ABGDS(J)  ! store parameters 
      END DO   
*     recursive loop for higher terms
      IAM=1
      DO K=2,NP
       IAM=3-IAM
       DO J=1,4
        ABGDS(J)=0.0D0
       END DO 
       DO I=1,N
        AL=AL+W(I)*X(I)*Q(3-IAM,I)**2       ! sum x*f(j-1)**2
        BE=BE+W(i)*X(I)*Q(IAM,I)*Q(3-IAM,I) ! sum f(j-2)*f(j-2)
       END DO 
       DO I=1,N
        Q(IAM,I)=(X(I)-AL)*Q(3-IAM,I)-BE*Q(IAM,I) ! unnormalized f(j) 
        GA=GA+W(I)*Q(IAM,I)*Q(IAM,I) ! sum for normalization
       END DO 
       IF(GA.NE.0.0) GA=1.0/DSQRT(GA)
       DO I=1,N
        Q(IAM,I)=GA*Q(IAM,I)   ! normalize f(j)
        DE=DE+W(I)*Q(IAM,I)*Y(I) ! sum for coefficient
       END DO 
       SQSUM=MAX(0.0D0,SQSUM-DE*DE) ! update sum of squares
       DO J=1,5
        P(J,K)=ABGDS(J)        ! store parameters
       END DO   
      END DO  
      P(1,1)=FLOAT(100*N+NP)   ! dimension parameter
      P(2,1)=1.0               ! factor for error calculation
      END

      SUBROUTINE ORTHFE(X,P,Y,DY)
*     Return function value Y and error DY of the orthogonal polynomial
*     at the abscissa value X, using the array P (see ORTHFT)                 
      REAL X,Y,DY,P(5,*)
      DOUBLE PRECISION Q(2),F,DF
      NP=MOD(P(1,1),100.0)+0.5
      Q(1)=P(3,1)       ! constant term
      Q(2)=0.0
      F =P(3,1)*P(4,1)  
      DF=P(3,1)*P(3,1)
      IAM=1
      DO K=2,NP
       IAM=3-IAM  ! recurrence relation for higher terms
       Q(IAM)=((X-P(1,K))*Q(3-IAM)-P(2,K)*Q(IAM))*P(3,K)
       F =F +P(4,K)*Q(IAM) 
       DF=DF+Q(IAM)*Q(IAM)
      END DO 
      Y =F               ! function value
      DY=SQRT(DF*P(2,1)) ! standard deviation 
      END