      program teach063
      implicit real*8 (a-h,o-z)
c
      n=1000
      pi=4*atan(1.0d0)
      g=9.8d0
      a=0.0d0
      T0=2*pi/sqrt(g)
      xmax=pi/2
      open(16,file='teach063.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 chuten(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 chuten(x0,a,b,n,sum)
c **********************************************************************
      implicit real*8 (a-h,o-z)
c
      sum = 0.0d0
      dx = (b - a)/n
      do i = 1, n
         x = a + (i - 0.5d0)*dx
         sum = sum + fv(x,x0)
      end do
      sum = sum*dx
C
      return
      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
