!###########################################################
! page.54 表3.4
! 台形公式，シンプソン公式，
! ガウス・ルジャンドル積分公式のプログラム
!
! S=(128.d0/3.d0) \int_0^1 dsqrt((x*(1.d0-x))**3) dx
!
! z=2*x-1 と変数変換すると
!
! x | 0 --> +1
! z |-1 --> +1
!
! dz=2.d0*dx
!
! であるから、
!
! S=(8.d0/3.d0) \int_{-1}^{+1} dsqrt(1.d0-z**2)**3 dz 
!
! を計算すれば良い。
!
!###########################################################
PROGRAM prog09
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: NMAX=10000
REAL(8),DIMENSION(0:NMAX) :: xx,yy
REAL(8),DIMENSION(NMAX) :: xx2,ww2
REAL(8),DIMENSION(0:NMAX) :: P
REAL(8) :: I1,I2,I3

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

WRITE(*,"(A)") 'N        台形公式          シンプソン公式    ガウス・ルジャンドル積分公式'

!********************
DO NS=1,10,1
IF(NS==1)THEN
N=2
ELSEIF(NS==2)THEN
N=3
ELSEIF(NS==3)THEN
N=4
ELSEIF(NS==4)THEN
N=5
ELSEIF(NS==5)THEN
N=6
ELSEIF(NS==6)THEN
N=10
ELSEIF(NS==7)THEN
N=20
ELSEIF(NS==8)THEN
N=100
ELSEIF(NS==9)THEN
N=1000
ELSEIF(NS==10)THEN
N=10000
ENDIF
!********************
xx=0.d0;yy=0.d0;xx2=0.d0;ww2=0.d0

h=2.d0/dble(N)
xx(0)=-1.d0
xx(N)=+1.d0
DO i=1,N-1
xx(i)=xx(0)+dble(i)*h
yy(i)=func(xx(i))
END DO

!**********************************************!
! 台形公式
!**********************************************!
I1=yy(0)+yy(N)
DO i=1,N-1
I1=I1+2.d0*yy(i)
END DO
I1=I1*h/2.d0

!**********************************************!
! シンプソン公式
!**********************************************!
M=N/2
I2=0.d0
DO i=0,M-1
I2=I2+yy(2*i)+4.d0*yy(2*i+1)+yy(2*i+2)
END DO
I2=I2*h/3.d0

!**********************************************!
! ガウス・ルジャンドル公式
!**********************************************!
CALL ZEROTEN(N,NMAX,xx2)

DO i=1,N
  a=0.d0
  CALL legendre(N,xx2(i),P,NMAX)
    DO k=0,N-1
    a=a+dble(2*k+1)*(P(k))**2
    END DO
  ww2(i)=2.d0/a
END DO

I3=0.d0
DO i=1,N
I3=I3+ww2(i)*func(xx2(i))
END DO


!********************
WRITE(*,"('N=',I5,3F18.14)") N,I1,I2,I3
!********************
END DO
!********************

STOP
END PROGRAM prog09

! ******************************************************************
!
!        ルジャンドル多項式
!       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 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!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION func(x)
   REAL(8),INTENT(IN) :: x
   REAL(8) :: func
   func=(8.d0/3.d0)*dsqrt((1.d0-x**2)**3)
RETURN
END FUNCTION func
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!