!###########################################################
! page.52 表3.1
! 台形公式，シンプソン公式，
! ガウス・ルジャンドル積分公式のプログラム
!
! f(x)=sqrt(1.d0-x**2)
!###########################################################
PROGRAM prog07
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: NMAX=200
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 N=8,20,4
!********************
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


!********************
E1=dabs(I1-pi/2.d0)
E2=dabs(I2-pi/2.d0)
E3=dabs(I3-pi/2.d0)
WRITE(*,"('N=',I3,1p3E16.8)") N,E1,E2,E3
!********************
END DO
!********************

STOP
END PROGRAM prog07

! ******************************************************************
!
!        ルジャンドル多項式
!       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=dsqrt(1.d0-x**2)
RETURN
END FUNCTION func
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!