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