!###########################################################
! page.35 表2.2
! チェビシェフ多項式補間のプログラム
!
! f(x)=1/(1+25*x**2)
!###########################################################
PROGRAM prog05
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: NMAX=50
REAL(8),DIMENSION(0:NMAX) :: ch
REAL(8),DIMENSION(0:NMAX) :: yy
REAL(8),DIMENSION(NMAX+1) :: bb
REAL(8),DIMENSION(NMAX+1,NMAX+1) :: aa

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

DO NS=1,3
WRITE(*,*)
!********************
DO N=12,48,4

yy=0.d0;bb=0.d0;aa=0.d0

h=2.d0/dble(N*ndx)

DO i=0,N
    x=funcx(i)
    yy(i)=func(x)
    bb(i+1)=yy(i)
	! チェビシェフ多項式
    CALL cheby(N,x,ch,NMAX)
      DO j=0,N
      aa(i+1,j+1)=ch(j)
      END DO
END DO

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



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

IF(NS==1)THEN
WRITE(*,'("Ez (",i2,")=",1PE15.8)') N+1,eps
ELSEIF(NS==2)THEN
WRITE(*,'("Ee1(",i2,")=",1PE15.8)') N+1,eps
ELSE
WRITE(*,'("Ee2(",i2,")=",1PE15.8)') N+1,eps
END IF

END DO
!********************
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

IF(NS==1)THEN
funcx=dcos(pi*dble(2*k+1)/dble(2*N+2))
ELSEIF(NS==2)THEN
funcx=dcos(pi*dble(k+1)/dble(N+2))
ELSE
funcx=dcos(pi*dble(k)/dble(N))
END IF

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 cheby(N,xx,ch,NMAX)
   funcy=0.d0
   DO k=0,N
   funcy=funcy+bb(k+1)*ch(k)
   END DO
RETURN
END FUNCTION funcy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

END PROGRAM prog05


! ******************************************************************
!
!        calculate chebyshev polynomials 
!
!      ch(i,j) ----- i-1 th dif. of the j th chebyshev polynomial
!
! ******************************************************************
SUBROUTINE cheby(N,x,ch,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) :: ch
ch=0.d0
ch(0)=1.d0
ch(1)=x
DO j=2,N
ch(j)=2.d0*x*ch(j-1)-ch(j-2)
END DO

RETURN
END SUBROUTINE cheby
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!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+1,NMAX+1) :: a
REAL(8),INTENT(INOUT),DIMENSION(NMAX+1) :: 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!