!*************************************************************
!*    Integration of a discrete function by the weighting    *
!*    coefficients Method                                    *
!* --------------------------------------------------------- *
!* SAMPLE RUN:                                               *
!*                                                           *
!* Number of points: 7                                       *
!*                                                           *
!* Input the 7 points, example: 1.5, 0.666666                *
!*                                                           *
!* 0  ?  1, 1                                                *
!* 1  ?  1.5, 0.666666                                       *
!* 2  ?  2, 0.5                                              *
!* 3  ?  2.5, 0.4                                            *
!* 4  ?  3, 0.333333                                         *
!* 5  ?  3.5, 0.285714                                       *
!* 6  ?  4, 0.25                                             *
!*                                                           *
!* Coefficients A(I):                                        *
!* A( 0 ) =  0.146428E+00                                    *
!* A( 1 ) =  0.771429E+00                                    *
!* A( 2 ) =  0.964289E-01                                    *
!* A( 3 ) =  0.971428E+00                                    *
!* A( 4 ) =  0.964289E-01                                    *
!* A( 5 ) =  0.771429E+00                                    *
!* A( 6 ) =  0.146429E+00                                    *
!*                                                           *
!* Integral of F(X) from  1.00 to   4.00                     *
!*   I = 0.138666E+01                                        *
!*                                                           *
!* --------------------------------------------------------- *
!* Reference: "Methodes de calcul numerique - Tome 1         *
!*             By Claude Nowakowski, PS1 1981" [BIBLI07].    *
!*                                                           *
!*                       F90 Release By J-P Moreau, Paris.   *
!*                              (www.jpmoreau.fr)            *
!*************************************************************
! See explanations in file dinteg.txt
! -----------------------------------
Program weighting_Coeffs

integer, parameter :: SIZE = 25 

real*8 X(0:SIZE), Y(0:SIZE), A(0:SIZE), T(0:SIZE, 0:SIZE), B(0:SIZE)

PRINT *,' '
WRITE(*,10,advance='no'); READ *, N

PRINT *,' '
WRITE(*,20)  N
PRINT *,' '

N = N - 1

DO K = 0, N
  WRITE(*,30,advance='no')  K
  READ *, X(K), Y(K)
END DO

! Calculate the coefficients of the linear system
DO J=0, N
  T(0, J) = 1.d0
  DO I = 1, N
    T(I, J) = T(I - 1, J) * X(J)
  END DO
END DO

! right side coefficients
AA = 1.d0; BB = 1.d0
DO I = 0, N
  AA = AA * X(0); BB = BB * X(N)
  B(I) = (BB - AA) / (I + 1)
END DO

! Calculate coefficients A(I) and triangularize
DO K = 0, N - 1
  DO I = K + 1, N
    B(I) = B(I) - T(I, K) / T(K, K) * B(K)
    DO J = K + 1, N
      T(I, J) = T(I, J) - T(I, K) / T(K, K) * T(K, J)
    END DO
  END DO
END DO

! Solve triangular linear system T * A = B
A(N) = B(N) / T(N, N)
DO I = N - 1, 0, -1
  S = 0.d0
  DO K = I + 1, N
    S = S + T(I, K) * A(K)
  END DO
  A(I) = (B(I) - S) / T(I, I)
END DO

PRINT *,' '
PRINT *,' Coefficients A(I):'
DO I = 0, N
  WRITE(*,40) I, A(I)
END DO

! Calculate Integral of discreet F(X)
S = 0.d0
DO I = 0, N
  S = S + A(I) * Y(I)
END DO

! Print final result
PRINT *,' '
WRITE(*,50) X(0), X(N)
WRITE(*,60) S
PRINT *,' '

10 format(' Number of points: ')
20 format(' Input the',I2,' points, example: 1.5, 0.666666')
30 format(' ', I2,'  ?  ')
40 format(' A(',I2,') = ', E12.6)
50 format(' Integral of F(X) from ', F5.2,' to ',F5.2)
60 format('   I = ',E12.6)

pause

END

!end of file dinteg.f90