      program teach064
      implicit real*8 (a-h,o-z)
c
      n=50
      pi=4*atan(1.0d0)
      g=9.8d0
      a=0.0d0
      T0=2*pi/sqrt(g)
      xmax=pi/2
      open(16,file='teach064.dat')
      write( *,'(a)') '# x0, T, T0, T-T0'
      write(16,'(a)') '# x0, T, T0, T-T0'
      do i=1,100
         x0=xmax/100*i
         b=x0
         call DEint(x0,a,b,n,DeltaT)
         T1=4.0d0/sqrt(g)*DeltaT+T0/sqrt(sin(x0))*sqrt(x0)
         write( *,800) real(x0),real(T1),real(T0),real(T1-T0)
         write(16,800) real(x0),real(T1),real(T0),real(T1-T0)
      enddo
 800  format(4(1x,1pe12.5))
      end
c **********************************************************************
      subroutine DEint(x0,a,b,n,sum)
c **********************************************************************
c
c     Integral of fv(x,x0) by using Double Exponetial Formula
c     by Takahashi-Mori.
c
c       sum= integral_{-infinity}^infinity FUNC(phi(t)) d(phi)/dt dt
c       phi = (b-a)/2*tanh(pi/2*sinh(t))+(b+a)/2
c
c **********************************************************************
      implicit real*8 (a-h,o-z)
c
      tmax=4.7d0
      dt=tmax/n
      pi=4.0d0*atan(1.0d0)
c
      sum = 0.0d0
      do i = -n,n
         t=i*dt
         s=sinh(t)
         phi=(b-a)/2*tanh(pi/2*s)+(b+a)/2
         dphidt=pi/2*cosh(t)/(cosh(pi/2*s))**2*(b-a)/2
         x=phi
         sum=sum+dphidt*fv(x,x0)
c        write(99,800) x,x0-x,cos(x)-cos(x0),sqrt(x0-x)
      end do
         sum = sum*dt
c800     format(4(1x,1pe19.12))
c        write(99,*)
 
      end
c **********************************************************************
      function fv(x,x0)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      eps=1.0d-12
      if(x.ge.x0) then
         fv=0.0d0
      elseif(x0-x.lt.eps) then
         fv=sqrt((x0-x)/2/sin(x0))/4*(1.0d0/tan(x0)-1.0d0/x0)
      else
         fv1=1.0d0/sqrt(2*(cos(x)-cos(x0)))
         fv0=1.0d0/sqrt(x0**2-x**2)/sqrt(sin(x0))*sqrt(x0)
         fv=fv1-fv0
      endif
      end
