!********************************************************
!*  LEAST SQUARE APPROXIMATION OF A DISCRETE FUNCTION   *
!*  USING ORTHOGONAL POLYNOMIALS                        *
!* ---------------------------------------------------- *
!* SAMPLE RUN:                                          *
!* (Approximate SIN(X) defined by 60 points, from x=0   *
!*  to x=PI/2, with degree K=7).                        *
!*                                                      *
!*  I     COEFFICIENTS   STD DEVIATION                  *
!*  0   0.63383395E+00  0.14142136E+00                  *
!*  1   0.66275983E+00  0.30570229E+00                  *
!*  2  -0.37764682E+00  0.73926175E+00                  *
!*  3  -0.11371820E+00  0.18223439E+01                  * 
!*  4   0.28616636E-01  0.45265886E+01                  *
!*  5   0.57493558E-02  0.11296326E+02                  *
!*  6  -0.96142085E-03  0.28296963E+02                  *
!*  7  -0.13769625E-03  0.71136688E+02                  *
!*                                                      *  
!*  I     VARIABLE        EXACT R         APPROX. R     *
!*  1   0.00000000E+00  0.00000000E+00 -0.26089835E-07  *
!*  2   0.31415927E+00  0.30901699E+00  0.30901698E+00  *
!*  3   0.62831853E+00  0.58778525E+00  0.58778525E+00  *
!*  4   0.94247780E+00  0.80901699E+00  0.80901700E+00  *
!*  5   0.12566371E+01  0.95105652E+00  0.95105651E+00  *
!*  6   0.15707963E+01  0.10000000E+01  0.99999997E+01  *
!*                                                      *
!* ---------------------------------------------------- *    
!* Ref.: From Numath Library By Tuan Dang Trong in      *
!*       Fortran 77 [BIBLI 18].                         *
!*                                                      *
!*              F90 Release 1.0 By J-P Moreau, Paris.   *
!*                        (www.jpmoreau.fr)             *
!********************************************************
      PARAMETER (K=7,M=50,NW=6)
      REAL*8 X(M),Y(M),SIGMA(M),S(0:K),ECART(0:K),ALPHA(K),  &
	  BETA(0:K-1),P,R,W,W0,DW,PI,WORK(4*M)
      PI=4.D0*ATAN(1.D0)
      DO I=1,M
      X(I)=PI*(I-1)/(2*(M-1))
      Y(I)=SIN(X(I))
      SIGMA(I)=1.D0
      ENDDO

      CALL  MCAPPR(K,M,X,Y,SIGMA,S,ALPHA,BETA,ECART,WORK)

	  PRINT *,' '
      WRITE(*,'(1X,A)') '  I     COEFFICIENTS    STD DEVIATION'
      WRITE(*,'(I4,2E17.8)') (I,S(I),ECART(I),I=0,K)
	  PRINT *,' '
      WRITE(*,'(1X,A)') '  I     VARIABLE         EXACT R          APPROX. R'
      W0=0.D0
      DW=PI/10.D0
      DO I=1,NW
      W=W0+(I-1)*DW
      R=P(K,S,ALPHA,BETA,W)
      WRITE(*,'(I4,3E17.8)') I,W,SIN(W),R
      ENDDO
      PRINT *,' '
      STOP
      END !OF MAIN PROGRAM

      SUBROUTINE MCAPPR(K,M,X,Y,SIGMA,S,ALPHA,BETA,ECART,WORK)
!=======================================================================
!     LEAST SQUARES APPROXIMATION OF A FUNCTION F(X) DEFINED BY M POINTS
!     X(I), Y(I) BY USING ORTHOGONAL POLYNOMIALS
!=======================================================================
!     INPUTS:
!     K   : DEGREE OF POLYNOMIALS
!     M   : NUMBER OF POINTS
!     X,Y : TABLES OF DIMENSION M TO STORE M ABSCISSAS AND
!           M ORDINATES OF GIVEN POINTS
!     SIGMA : TABLE OF DIMENSION M TO STORE THE STANDARD DEVIATIONS
!             OF VARIABLE Y
!     OUTPUTS:
!     S    : TABLE OF DIMENSION(0:K)
!     ALPHA: TABLE OF DIMENSION (K)
!     BETA : TABLE OF DIMENSION (0:K-1)
!     WORKING SPACE:
!     WORK : TABLE OF DIMENSION (4*M)
!     NOTE:
!     COEFFICIENTS S,ALPHA,BETA ARE USED TO EVALUATE VALUE
!     AT POINT Z OF THE BEST POLYNOMIAL OF DEGREE K
!     BY USING FUNCTION P(K,S,ALPHA,BETA,Z) DESCRIBED BELOW.
!=======================================================================
      REAL*8 X(M),Y(M),SIGMA(M),S(0:K),ECART(0:K),ALPHA(K),  &
	  BETA(0:K-1), WORK(4*M)
      M1=M+1
      M2=M+M1
      M3=M+M2
      CALL MCARRE(K,M,X,Y,SIGMA,S,ALPHA,BETA,ECART,WORK(1),WORK(M1),  &
                  WORK(M2),WORK(M3))
      END

      !LEAST SQUARES SUBROUTINE
      SUBROUTINE MCARRE(K,M,X,Y,SIGMA,S,ALPHA,BETA,ECART,P1,P2,P3,P4)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(M),Y(M),SIGMA(M),S(0:K),ECART(0:K),ALPHA(K),  &
      BETA(0:K-1),P1(M),P2(M),P3(M),P4(M)
      DO I=1,M
      P1(I)=0.D0
      P2(I)=1.D0
      ENDDO
      W=M
      BETA(0)=0.D0
      DO I=0,K-1
      OMEGA=PRD(Y,P2,SIGMA,M)
      S(I)=OMEGA/W
      T=PRD(P2,P2,SIGMA,M)/W**2
      ECART(I)=SQRT(T)
      DO L=1,M
      P4(L)=X(L)*P2(L)
      ENDDO
      ALPHA(I+1)=PRD(P4,P2,SIGMA,M)/W
      DO L=1,M
      P3(L)=(X(L)-ALPHA(I+1))*P2(L)-BETA(I)*P1(L)
      ENDDO
      WPR=PRD(P3,P3,SIGMA,M)
      IF(I+1.LE.K-1) BETA(I+1)=WPR/W
      W=WPR
      DO L=1,M
      P1(L)=P2(L)
      P2(L)=P3(L)
      ENDDO
      ENDDO
      OMEGA=PRD(Y,P2,SIGMA,M)
      S(K)=OMEGA/W
      T=PRD(P2,P2,SIGMA,M)/W**2
      ECART(K)=SQRT(T)
      RETURN
      END
      FUNCTION PRD(X,Y,Z,M)
      REAL*8 X(M),Y(M),Z(M),PRD,SUM
      SUM=0.D0
      DO I=1,M
      SUM=SUM+X(I)*Y(I)/Z(I)**2
      ENDDO
      PRD=SUM
      RETURN
      END

      FUNCTION P(K,S,ALPHA,BETA,X)
!=====================================================================
!     THIS FUNCTIONE ALLOWS EVALUATING VALUE AT POINT X OF A FUNCTION
!     F(X), APPROXIMATED BY A SYSTEM OF ORTHOGONAL POLYNOMIALS Pj(X)
!     THE COEFFICIENTS OF WHICH, ALPHA,BETA HAVE BEEN DETERMINED BY
!     LEAST SQUARES.
!=====================================================================
      REAL*8 S(0:*),ALPHA(*),BETA(0:*),P,B,X,BPR,BSD
      B=S(K)
      BPR=S(K-1)+(X-ALPHA(K))*S(K)
      DO I=K-2,0,-1
      BSD=S(I)+(X-ALPHA(I+1))*BPR-BETA(I+1)*B
      B=BPR
      BPR=BSD
      ENDDO
      P=BPR
      RETURN
      END

!end of file approx1.f90