!###########################################################
! page.56 表3.5
! 台形公式のプログラム
! f(x)=exp(-x**2)
!
! I=\int_{-L}^{L} f(x) dx
!
!###########################################################
PROGRAM prog10
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: NMAX=10000
REAL(8),DIMENSION(0:NMAX) :: xx,yy
REAL(8) :: I1,LL
REAL(8),DIMENSION(10) :: E3
CHARACTER(3) :: z='   '

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

WRITE(*,"(A,9I6)") 'L=   ',(i,i=10,2,-1)
WRITE(*,'(100A)') ('-',i=1,60)
!********************
DO NH=1,14,1
h=dble(NH)/10.d0
E3=0.d0
DO L=10,2,-1
!********************
xx=0.d0;yy=0.d0

! hとLからNを計算する
! NINT は最も近い整数を返す関数
! ここでは，hとNからL(整数)を計算するとき，
! L(整数)に最も近くなるようなN(整数)を計算している。

N=NINT(2.d0*dble(L)/h)

! N と h から LL(実数)を計算し，
! -LL〜LLを積分範囲とする
! L(整数)とは若干異なる
LL=dble(N)*h/2.d0
xx(0)=-LL
xx(N)=LL

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

E1=dabs(I1-dsqrt(pi))

IF(E1==0.d0)THEN
E3(L)=0.d0
ELSE
  DO ie=0,15
  IF((E1<10.d0**(-ie)).and.(E1>=10.d0**(-(ie+1))))THEN
  E2=E1*10.d0**(ie+1)
  E3(L)=E1/E2
  EXIT
  END IF
  END DO
ENDIF

END DO
WRITE(*,"('h=',F3.1,9(1PE6.0,TL6,A,TR3))") h,(E3(i),z,i=10,2,-1)
END DO
!********************

WRITE(*,*)
WRITE(*,"(A)") '注意 : +00 の部分は誤差が 0 (本文中では -K の部分)'

STOP
END PROGRAM prog10

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION func(x)
   REAL(8),INTENT(IN) :: x
   REAL(8) :: func
   func=dexp(-x**2)
RETURN
END FUNCTION func
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!