!*****************************************************************************!
! page53 表3.2
! ガウス・ルジャンドル積分公式のプログラム
! f(x)=1.d0//sqrt(1.d0-x**2)
! x=[-1,1]
!*****************************************************************************!

!*****************************************************************************!
! パラメータ
!*****************************************************************************!
module param
   real(8) :: pi,eps,s0
end module param
!*****************************************************************************!
!  メインプログラム
!*****************************************************************************!
program main
   use param
   implicit none
   integer N
   real(8) x(1:200),w(1:200),y(1:200)
   real(8) s

   call para()
   do N=50,200,50
      call zeroten(N,x)
      call weight(N,w,x)
      call sekibun_calc(N,x,y,w,s)  
      write(*,*) '|I-π|', dabs(s0-s)
   end do

   stop
end program main
!*****************************************************************************!
! パラメータの設定
! eps : 二分法の収束判定
! s0  : 求積解
!*****************************************************************************!
subroutine para()
   use param
   implicit none

   pi=4.d0*datan(1.d0)
   eps=1.d-15
   s0=pi

   return
end subroutine para
!*****************************************************************************!
!  ゼロ点を求める(二分法)
!*****************************************************************************!
subroutine zeroten(nn,x)
   use param
   implicit none
   integer i,m,nn
   real(8) a,b,c,fa,fb,fc
   real(8) x(1:200)

   m=int( dble(nn/2.0d0) )
   x(m)=0.0d0  ! Nが奇数のためのもの
   do i=1, m
      a=dsin( pi*dble(nn-2*i  )/dble(2*nn) )
      b=dsin( pi*dble(nn-2*i+2)/dble(2*nn) )
      call pnx(nn, a, fa)
      call pnx(nn, b, fb)
      if(fa*fb > 0.0d0) then
         write(*,*) 'Error a,b, i=' , i
         exit
      end if
      fc=1.0d0
      do while( dabs(fc) > eps .and. dabs(a-b) > eps)
         c=(a+b)/2.0d0
         call pnx(nn, c, fc)
         if(fc*fa < 0.0d0) then
            b=c
         else
            a=c
         end if
      end do
      x(nn-i+1)=c
      x(i)=-c
   end do
	
   return
end subroutine zeroten
!*****************************************************************************!
!  ルジャンドル多項式
!*****************************************************************************!
subroutine pnx(nn,abc,fabc)
   use param
   implicit none
   integer j,nn
   real(8) p0,p1,p2,abc,fabc

   p0=1.0d0
   p1=abc
   if(nn==0) then
      fabc=p0
      return
   end if
   if(nn==1) then
      fabc=p1
      return
   end if
   do j=1, nn-1
      p2=( dble(2*j+1)*abc*p1 - dble(j)*p0 )/dble(j+1)
      p0=p1
      p1=p2
   end do
   fabc=p2

   return
end subroutine pnx
!*****************************************************************************!
! 重みを計算
!*****************************************************************************!
subroutine weight(nn,w,x)
   use param
   implicit none
   integer i,k,nn
   real(8) f,pnk
   real(8) x(1:200),w(1:200)

   do i=1, nn
      f=0.0d0
      do k=0, nn-1
         call pnx(k, x(i), pnk)
         f=f+(2.0d0*dble(k)+1.0d0)*pnk*pnk
      end do
      w(i)=2.0d0/f
   end do

   return
end subroutine weight
!*****************************************************************************!
! 与えられた関数yを積分
!*****************************************************************************!
subroutine sekibun_calc(nn,x,y,w,integral)
   use param
   implicit none 
   integer i, nn
   real(8) y(1:200),w(1:200),x(1:200)
   real(8) integral

   do i=1, nn
      y(i)=1.d0/dsqrt(1.0d0-x(i)*x(i))
   end do

   integral=0.0d0
   do i=1, nn
      integral=integral+w(i)*y(i)
   end do

   return
end subroutine sekibun_calc