!###########################################################
! page.60 表3.7
! 台形公式のプログラム(座標変換)
!
! f(x)=1.d0/dsqrt(1.d0-x**2)
! I=\int_{-1}^{1} f(x) dx
!
! を、x=tanh(z) と変数変換すると
!
! I=\int_{-\infty}^{\infty} 1.d0/cosh(z) dz
!
!###########################################################
PROGRAM prog12
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(20) :: E3
INTEGER,DIMENSION(20) :: N3
CHARACTER(3) :: z='   '

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

WRITE(*,"(A,9I11)") 'L=   ',(i,i=20,15,-5)
WRITE(*,'(100A)') ('-',i=1,30)
!********************
DO NH=1,14,1
h=dble(NH)/10.d0
E3=0.d0
N3=0.d0
DO L=20,15,-5
!********************
xx=0.d0;yy=0.d0

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

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

! 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-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,'(',I3,')'))") h,(E3(i),z,N3(i),i=20,15,-5)
END DO
!********************

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

STOP
END PROGRAM prog12

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION func(x)
   REAL(8),INTENT(IN) :: x
   REAL(8) :: func
   func=1.d0/dcosh(x)
RETURN
END FUNCTION func
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!