program prog033
      implicit real*8 (a-h,o-z)
c
c     This program gives the integral
c     integral_0^1 FUNC(x) dx
c     by using the Double Exponential Formula
c     by Takahashi-Mori
c
c     Functions integrated in this program are
c
c        f0(x)=1.0d0/sqrt(x)
c        f1(x)=pi*sin(pi*x)
c        f2(x)=4*sqrt(1.0d0-x**2)
c        f3(x)=-log(x)
c
      external f1,f2,f3,f0
      a=0.0d0
      b=1.0d0
c
      hmax=4.7d0
c
      exa1=2.0d0
      exa2=4.0d0*atan(1.0d0)
      exa3=1.0d0
c
      write( *,'(a)') '#   n  dx          Err1        Err2        Err3'
      write(16,'(a)') '#   n  dx          Err1        Err2        Err3'
      do n=5,50,5
         call DEint(f1,a,b,n,sum1)
         call DEint(f2,a,b,n,sum2)
         call DEint(f3,a,b,n,sum3)
c
         err1=abs(sum1-exa1)
         err2=abs(sum2-exa2)
         err3=abs(sum3-exa3)
c
      write( *,800) n,hmax/n,err1,err2,err3
      write(16,800) n,hmax/n,err1,err2,err3
 800  format(1x,i4,4(1x,g11.4))
      enddo
c
      end
c **********************************************************************
      subroutine DEint(FUNC,a,b,n,sum)
c **********************************************************************
c
c     This subroutine calculate the integral
c       sum = integral_a^b FUNC(x) dx
c     by using the Double Exponential formula (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)
      external FUNC
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
         sum=sum+FUNC(phi)*dphidt
      end do
         sum = sum*dt
c
      end
c **********************************************************************
      function f0(x)
c **********************************************************************
c     function f0(x)=1.0d0/sqrt(x)
c        0 < x < eps it returns 1.0d0/sqrt(eps)
c        x < 0       it returns 0
c **********************************************************************
      implicit real*8 (a-h,o-z)
c
      eps=1.0d-32
      if(x.lt.0) then
         f0=0.0d0
      else if(x.lt.eps) then
         f0=1.0d0/sqrt(eps)
      else
         f0=1.0d0/sqrt(x)
      endif
c
      end
c **********************************************************************
      function f1(x)
c **********************************************************************
c     function f1(x)=pi*sin(pi*x)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      pi=4.0d0*atan(1.0d0)
c
      f1=pi*sin(pi*x)
c
      end
c **********************************************************************
      function f2(x)
c **********************************************************************
c     function f2(x)=4*sqrt(1.0d0-x**2)
c        |x| > 1          it returns 0
c **********************************************************************
      implicit real*8 (a-h,o-z)
c
      if(abs(x).gt.1.0d0) then
         f2=0.0d0
      else
         f2=4*sqrt(1.0d0-x**2)
      endif
c
      end
c **********************************************************************
      function f3(x)
c **********************************************************************
c     function f3(x)=-log(x)
c        0 < x < eps  it returns -log(eps)
c        x < 0        it returns 0
c **********************************************************************
      implicit real*8 (a-h,o-z)
c
      eps=1.0d-32
      if(x.lt.0) then
         f3=0.0d0
      else if(x.lt.eps) then
         f3=-log(eps)
      else
         f3=-log(x)
      endif
c
      end