!***************************************************
!*  Evaluate a Legendre Polynomial P(x) of Order n *
!*  for argument x by using Horner's Rule          *
!* ----------------------------------------------- *
!* SAMPLE RUN:                                     *
!*                                                 *
!*  give order: 7                                  *
!*                                                 *
!*  give argument x: 0.5                           *
!*                                                 *
!*   P(x) =   0.2231445                            *
!*                                                 *
!* ----------------------------------------------- *
!*               F90 Release By J-P Moreau, Paris. *
!*                       (www.jpmoreau.fr)         *
!***************************************************
Program Eval_Leg

Implicit None

integer, parameter :: NMAX = 12

real A(0:NMAX)
integer n
real x, Eval

  print *,' '
  write(*,10,advance='no'); read *, n
  print *,' '
  write(*,20,advance='no'); read *, x

  Call Legendre_Coeff(n, A)

  print *,' '
  print *,' P(x) = ', Eval(n,A,x)
  print *,' '

10 format(' give order: ')
20 format(' give argument x: ') 

End


!******************************************************
!* Legendre series coefficients evaluation subroutine *
!* by means of recursion relation. The order of the   *
!* polynomial is n. The coefficients are returned in  *
!* A(i) by increasing powers.                         *
!******************************************************
Subroutine Legendre_Coeff(n, A)
  Implicit None
  integer, parameter :: NMAX = 12
  integer n, i,j
  real A(0:NMAX)
  real B(0:NMAX,0:NMAX)
! Establish p0 and p1 coefficients
  B(0,0)=1.0; B(1,0)=0.0; B(1,1)=1.0
! Return if order is less then 2
  if (n > 1) then
	do i=2, n
      B(i,0) = -Real(i-1)*B(i-2,0)/Real(i)
      do j=1, i
        !Basic recursion relation
        B(i,j)=Real(i+i-1)*B(i-1,j-1)-Real(i-1)*B(i-2,j)
        B(i,j) = B(i,j) / Real(i)
      end do
    end do
    do i=0, n
	  A(i)=B(n,i)
    end do
  end if
  return
End


Subroutine HORNER (N, A, X0, K, B, R)
!--------------------------------------------------------------------
!     EVALUATE A POLYNOMIAL AND ITS DERIVATIVES BY HORNER'S METHOD
!
!     INPUTS:
!     N       ORDER OF POLYNOMIAL
!     A       VECTOR OF SIZE N+1 STORING THE COEFFICIENTS OF
!             POLYNOMIAL IN DECREASING ORDERS OF POWERS
!             I.E. P(X) = A(1)*X**N+A(2)*X**(N-1)+...+A(N)*X+A(N+1)
!     X0      GIVEN ARGUMENT
!     K       MAXIMUM ORDER OF DERIVATIVES ASKED FOR
!     B       FLAG
!             = .TRUE.  EVALUATION ONLY
!             = .FALSE. DETERMINATION OF POLYNOMIAL COEFFICIENTS
!                NEAR X0
!     R       VECTOR OF SIZE K+1 CONTAINS:
!             -  VALUES P(X0), P'(X0), P''(X0),..PK(X0),
!                IF B IS TRUE AND IF K < N
!             -  THE N+1 COEFFICIENTS OF POLYNOMIAL
!                Q(X-X0) = R(1)+R(2)*(X-X0)+...+R(N+1)*(X-X0)**(N),
!                IF B IS FALSE AND IF K = N.
!
!     REFERENCE:
!     ALGORITHM 337, COLLECTED ALGORITHMS FROM CACM, W.PANKIEWICKZ
!---------------------------------------------------------------------
  Implicit None
  integer N, K
  real A(*), R(*), X0
  logical B
  real RR
  integer I,J,L,NMJ

  RR = A(1)
  do I=1, K+1
    R(I) = RR
  end do
  do J=2, N+1
    R(1) = R(1)*X0 + A(J)
    NMJ = N-J+1
    if (NMJ > K) then
      L = K
    else
      L = NMJ
    end if
    do I=2, L+1
	  R(I) = R(I)*X0 + R(I-1)
    end do
  end do
  if (B==.TRUE.) then
    L = 1
    do I=2, K+1
      L = I*(I-1)
      R(I) = R(I)*L
    end do
  end if
  Return
End

! Evaluate P(x) of order n for argument x
real Function Eval(n, A, x)
  Implicit None
  integer n
  real A(0:n), x
  integer i
  integer, parameter :: NMAX = 12
  real A1(1:NMAX), R(1:NMAX)
  do i=n, 0, -1
    A1(n-i+1) = A(i)
  end do
  Call HORNER(n,A1,x,0,.TRUE.,R)
  Eval = R(1)
  Return
End

! end of file eval_leg.f90