!###########################################################
! page.40 図2.10
! ルジャンドル多項式補間のプログラム
!
! f(x)=1/(1+25*x**2)
!###########################################################
PROGRAM prog06
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: NMAX=50
REAL(8),DIMENSION(0:NMAX) :: P
REAL(8),DIMENSION(NMAX) :: xx,yy
REAL(8),DIMENSION(NMAX) :: bb
REAL(8),DIMENSION(NMAX,NMAX) :: aa

pi=4.d0*datan(1.d0)
ndx=100

!********************
DO NS=1,3
IF(NS==1)THEN
N=5
ELSEIF(NS==2)THEN
N=9
ELSE
N=17
ENDIF
IF(NS==1)THEN
OPEN(1,file='lege01.dat',status='unknown')
OPEN(2,file='lege00.dat',status='unknown')
ELSEIF(NS==2)THEN
OPEN(1,file='lege02.dat',status='unknown')
ELSE
OPEN(1,file='lege03.dat',status='unknown')
ENDIF

xx=0.d0;yy=0.d0;bb=0.d0;aa=0.d0
h=2.d0/dble(N*ndx)


CALL ZEROTEN(N,NMAX,xx)

DO i=1,N
    yy(i)=func(xx(i))
    bb(i)=yy(i)
	! ルジャンドル多項式
    CALL legendre(N,xx(i),P,NMAX)
      DO j=0,N-1
      aa(i,j+1)=P(j)
      END DO
END DO

! ガウスの消去法
CALL gauss(N,aa,bb,NMAX)

eps=0.d0
DO i=0,N-1
DO ii=1,ndx
x=-1.d0+dble(i*ndx+ii)*h
ylege=funcy(x)
ytrue=func(x)
eps=eps+(ytrue-ylege)**2

      WRITE(1,*) x,ylege
      IF(NS==1)THEN
      WRITE(2,*) x,ytrue
      ENDIF
END DO
END DO
eps=eps*h

WRITE(*,'("E(N=",i2,")=",1PE15.8)') N,eps

CLOSE(1,status='keep')
CLOSE(2,status='keep')
END DO
!********************

STOP
!###########################################################
CONTAINS

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION func(xx)
   REAL(8),INTENT(IN) :: xx
   REAL(8) :: func
   func=1.d0/(1.d0+25.d0*xx**2)
RETURN
END FUNCTION func
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION funcx(k)
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)
   INTEGER,INTENT(IN) :: k
   REAL(8) :: funcx

funcx=dcos(pi*dble(2*k+1)/dble(2*N+2))

RETURN
END FUNCTION funcx
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION funcy(xx)
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)
   REAL(8),INTENT(IN) :: xx
   REAL(8) :: funcy
call legendre(N,xx,P,NMAX)
   funcy=0.d0
   DO k=0,N
   funcy=funcy+bb(k+1)*P(k)
   END DO
RETURN
END FUNCTION funcy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

END PROGRAM prog06


! ******************************************************************
!
!        ルジャンドル多項式
!       P_j(x) (j=0,1,2,...,N)
!
! ******************************************************************
SUBROUTINE legendre(N,x,P,NMAX)
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)
INTEGER,INTENT(IN) :: N,NMAX
REAL(8),INTENT(IN) :: x
REAL(8),INTENT(OUT),DIMENSION(0:NMAX) :: P
P=0.d0
P(0)=1.d0
P(1)=x
DO j=1,N-1
P(j+1)=(dble(2*j+1)*x*P(j)-dble(j)*P(j-1))/dble(j+1)
END DO
RETURN
END SUBROUTINE legendre
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo!
SUBROUTINE gauss(n,a,b,NMAX)
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)
INTEGER,INTENT(IN) :: n,NMAX
REAL(8),INTENT(INOUT),DIMENSION(NMAX,NMAX) :: a
REAL(8),INTENT(INOUT),DIMENSION(NMAX) :: b
REAL(8),DIMENSION(n) :: ip,wk
eps=3.52d-15
ip=0.d0;wk=0.d0

DO k=1,n-1
!  find maximum element in the k-th column.
amax=dabs(a(k,k))
ipk=k

  DO i=k+1,n
  aik=dabs(a(i,k))
    IF(aik>amax) THEN
    ipk=i
    amax=aik
    END IF
  END DO

ip(k)=ipk

!  matrix is singular.
IF(amax<eps) THEN
WRITE(*,*)'(subr.glu) matrix is singular at k =',k
STOP
END IF

      IF(ipk/=k) THEN
      w=a(ipk,k)
      a(ipk,k)=a(k,k)
      a(k,k)=w
      END IF
!  compute alfa
      DO i=k+1,n
      a(i,k)=-a(i,k)/a(k,k)
      wk(i)=a(i,k)
      END DO
      
      DO j=k+1,n
         IF(ipk/=k) THEN
         w=a(ipk,j)
         a(ipk,j)=a(k,j)
         a(k,j)=w
         END IF
!  gaussian elimination
      t=a(k,j)
         DO i=k+1,n
         a(i,j)=a(i,j)+wk(i)*t
         END DO
      END DO
END DO
!  right-hand side
!  forward elimination process
DO k=1,n-1
   IF(ip(k)/=k) THEN
     w=b(ip(k))
     b(ip(k))=b(k)
     b(k)=w
   END IF

   t=b(k)

   DO i=k+1,n
   b(i)=b(i)+a(i,k)*t
   END DO
END DO

!  backward substitution process
DO k=n,2,-1
   b(k)=b(k)/a(k,k)
      t=b(k)
      DO i=1,k-1
      b(i)=b(i)-a(i,k)*t
      END DO
END DO
b(1)=b(1)/a(1,1)

RETURN
END SUBROUTINE gauss
!oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo!
!oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo!
SUBROUTINE ZEROTEN(N,NMAX,xx)
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)
INTEGER,INTENT(IN) :: N,NMAX
REAL(8),INTENT(OUT),DIMENSION(NMAX) :: xx
REAL(8),DIMENSION(0:NMAX) :: P

pi=4.d0*datan(1.d0)
xx=0.d0

delta=1.d-10
eps=1.d-10
M=N/2

DO i=1,M

a=dsin(pi*dble(N-2*i)/dble(2*N))
b=dsin(pi*dble(N-2*i+2)/dble(2*N))

CALL legendre(N,a,P,NMAX)
fa=P(N)
CALL legendre(N,b,P,NMAX)
fb=P(N)

IF(fa*fb>0) THEN
WRITE(*,*) 'Error a,b,  i==',i
ENDIF

fc=1.d0

DO WHILE(dabs(fc)>delta.and.dabs(a-b)>eps)
    diff=dabs(a-b)
    c=(a+b)/2.d0
    CALL legendre(N,c,P,NMAX)
    fc=P(N)

    IF(fc*fa>0.d0)THEN
    a=c ; fa=fc
    ELSE
    b=c ; fb=fc
    END IF
END DO

100 CONTINUE
xx(N-i+1)=c
xx(i)=-c

END DO

RETURN
END SUBROUTINE ZEROTEN
!oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo!